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Director’s  Foreword 


The  U.S.  Army  Research  Laboratory  (ARL)  mission  is  to  “Provide  innovative  science,  technology,  and 
analyses  to  enable  full  spectrum  operations.”  As  the  Army’s  corporate  laboratory,  we  provide  the 
technological  underpinnings  critical  to  providing  capabilities  required  by  our  current  and  future 
Soldiers. 

Our  nation  is  projected  to  experience  a  shortage  of  scientists  and  engineers.  ARL  recognizes  the 
criticality  of  intellectual  capital  in  generating  capabilities  for  the  Army.  As  the  Army’s  corporate 
laboratory,  addressing  the  projected  shortfall  is  a  key  responsibility  for  us.  We  have,  therefore, 
identified  the  nation’s  next  generation  of  scientists  and  engineers  as  a  key  community  of  interest  and 
have  generated  a  robust  educational  outreach  program  to  strengthen  and  support  them.  We  have 
achieved  many  successes  with  this  community.  We  believe  that  the  breadth  and  depth  of  our  outreach 
programs  will  have  a  significant  positive  effect  on  the  participants,  facilitating  their  journey  toward 
becoming  this  Nation’s  next  generation  of  scientists  and  engineers. 

A  fundamental  component  of  our  outreach  program  is  to  provide  students  research  experiences  at  ARL. 
During  the  summer  of  2010,  we  supported  research  experiences  at  ARL  for  over  150  undergraduate  and 
graduate  students.  Each  of  these  students  writes  a  paper  describing  the  results  of  the  work  they 
performed  while  at  ARL.  All  of  the  papers  were  of  high  quality,  but  only  a  few  could  be  presented  at 
our  student  symposium.  The  abstracts  for  all  papers  prepared  this  summer  are  contained  in  this  volume 
of  the  proceedings  and  they  indicate  that  there  were  many  excellent  research  projects  with  outstanding 
results.  It  is  unfortunate  that  there  was  not  enough  time  for  us  to  have  all  of  the  papers  presented.  We 
would  have  enjoyed  hearing  them  all. 

We  are  very  pleased  to  have  hosted  this  outstanding  group  of  students  for  the  summer.  It  is  our  hope 
that  they  will  continue  their  pursuit  of  technical  degrees  and  will  someday  assist  us  in  providing  critical 
technologies  for  our  Soldiers. 
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Introduction 


The  ARL  Summer  Student  Research  Symposium  is  an  ARL  Director’s  Award  Program  for  all  the 
students  participating  in  various  summer  scholarship  and  contract  activities  across  ARL.  The  goal  of  the 
program  is  to  recognize  and  publicize  exceptional  achievements  made  by  the  students  and  their  mentors 
in  the  support  of  Army  science. 

All  college  undergraduate  and  graduate  students  receiving  research  appointments  and  conducting  summer 
studies  at  ARL  are  automatically  enrolled  in  the  symposium  program.  As  an  integral  part  of  their  summer 
study,  all  students  are  expected  to  write  a  paper  on  their  work  which  summarizes  their  major  activity  and 
its  end  product. 

The  program  is  conducted  on  two  separate  competitive  levels:  undergraduate  and  graduate.  The  format 
of  the  paper  in  both  levels  is  the  same.  However,  the  evaluation  will  take  into  consideration  the 
difference  in  the  academic  level  of  the  students. 

All  students  submitted  their  research  paper  for  directorate  review.  Directorate  judging  panels  selected 
one  or  two  papers  from  each  competition  category  for  the  laboratory-wide  competition  at  the  Summer 
Student  Symposium  on  10  August  2010. 

Students  selected  by  their  directorate  for  competition  participated  in  the  one-day  Summer  Student 
Symposium  on  10  August  2010.  At  the  symposium,  the  students  presented  their  papers  to  an  audience  of 
ARL  scientists  and  engineers,  including  the  ARL  Director  and  an  ARL  Fellows  panel. 

This  volume  of  the  Summer  Student  Symposium  Proceedings  contains  the  papers  presented  at  the 
symposium. 
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Abstract 


The  concentration  of  energy  and  creation  of  hot  spots  in  crystalline  hexahydro-l,3,5-trinitro- 
1,3,5-s-triazine  (RDX)  is  influenced  by  localized  deformation  features.  This  work  presents 
results  from  molecular  dynamics  simulations  of  RDX  using  a  flexible  molecule  potential  energy 
function  to  elucidate  the  formation  of  these  localized  features.  Equilibration  simulations  are  used 
to  validate  the  orthotropic  a-  and  y-RDX  crystal  structure  and  properties  to  experiment.  The 
competition  between  cleavage  fracture  and  dislocation  nucleation  at  a  crack  tip  in  a-RDX  is  then 
evaluated  using  Rice’s  continuum  formulation.  Rice’s  analysis  leads  to  a  new  material 
parameter,  the  unstable  stacking  fault  energy,  which  is  the  maximum  energy  encountered  for 
sliding  on  a  slip  plane  from  molecular  dynamics  simulations.  This  work  determines  the  gamma 
surface,  or  stacking  energy,  on  the  (010)  and  (001)  slip  planes.  Maximum  saddle  points  from  the 
gamma  surface  are  compared  to  the  free  surface  energy  to  determine  the  likely  nature  of 
deformation  on  these  planes. 
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1.  Introduction/Background 


The  first  step  towards  developing  insensitive  munitions  is  uncovering  and  understanding  the 
subcontinuum  features  and  mechanisms  of  the  constituent  materials  that  make  up  the  energetic 
fill.  Identifying  structure-property  relationships  and  mechanisms  associated  with  inelastic 
response  due  to  loads  is  critical  for  building  up  a  first-principles  description  of  initiation  and 
failure,  which,  in  turn,  can  help  designers  and  engineers  develop  munitions  that  are  more 
effective  in  their  missions  but  safer  to  handle.  Today,  owing  to  the  confluence  of  high  rate 
physics,  interdisciplinary  effects,  and  the  imperfections  of  the  material  at  atomistic  scales,  it  still 
remains  unknown  precisely  what  leads  to  initiation  of  energetic  materials.  Though  many  events, 
both  deterministic  and  stochastic,  are  believed  to  be  involved,  their  precise  sequence  and  relative 
contributions  to  the  final  outcome  are  still  largely  unknown. 

The  energetic  fill  is  a  composite  structure  of  energetic  crystals,  such  as  crystalline  hexahydro- 
l,3,5-trinitro-l,3,5-s-triazine  (RDX),  encased  in  a  plastic  binder  material.  It  is  generally  agreed 
upon  that  decomposition  of  the  energetic  fill  is  the  culmination  of  micron-sized  reactions,  called 
hot  spots,  that  result  from  shockwave  interactions  and  dynamic  localization  phenomena  in  the 
heterogeneous  composition  of  the  explosive  fill.  Heterogeneities  in  the  fill  and  its  composite 
structure  come  in  the  form  of  voids,  grain  boundaries,  multimaterial  interfaces,  and  crystal 
orientation  mismatches.  The  mechanism  leading  to  localization  of  energy  and  decomposition  in 
single  crystals  of  energetic  material  free  of  heterogeneities  is  unclear.  Armstrong  (3)  proposes 
that  pile-ups  of  dislocations  are  suddenly  released  as  an  avalanche  of  localized  deformation  and 
available  energy  for  adiabatic  heating,  leading  to  the  nucleation  of  hot  spots.  For  the  energetic 
crystal  pentaerythritol  tetranitrate  (PETN),  Dick  and  Ritchie  ( 4 )  suggested  that  the  complex 
structure  of  the  molecule  sterically  hindered  displacement  on  slip  planes  and  led  to  deformed 
molecules,  whose  bonds  were  easily  broken.  Both  models  are  not  mutually  exclusive:  steric 
hindrance  may  aid  in  the  pile  up  of  dislocations  where  the  result  could  be  dislocation  avalanches 
and  fracture,  or  planes  of  highly  resolved  shear  stresses  and  deformed  molecules.  However, 
there  is  a  lack  of  experimental  observations  made  of  the  deformation  mechanisms  in  molecular 
crystals  to  provide  conclusive  evidence  as  to  which  mechanism  is  dominant. 

In  this  work,  molecular  dynamics  will  be  used  to  study  the  atomistic  deformation  features  of 
RDX  to  provide  a  mechanistic  basis  for  its  brittle  nature  and  a  description  of  operable 
deformation  mechanisms.  We  will  present  calculations  and  their  analyses  that  clarify  the 
competitive  balance  between  dislocation  nucleation  and  cleavage  at  crack  tips  in  RDX. 
Dislocation  emission  competes  with  Griffith  cleavage  near  crack  tips  in  much  the  same  way  as 
ductility  and  brittleness  in  metals  have  long  been  studied.  Limited  dislocation  activity  at  the 
crack  tip  becomes  entangled  by  the  low  symmetry  crystal,  and  complex  molecule  shape  and 
subsequent  deformation  may  lead  to  an  avalanche  of  localized  dislocations  and  fracture  or  large 
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shear  stresses  sufficient  for  molecule  conformational  changes.  Using  nanoindentation,  Ramos  et 
al.  (5)  made  experimental  observations  of  limited  plastic  slip  on  the  slip  plane,  followed  by 
fracture  at  low  loads.  From  their  work,  they  were  able  to  determine  the  slip  systems,  but  the 
cause  of  the  low  fracture  threshold  was  unresolved.  Atomistic  level  simulations  of  these 
experimental  slip  and  cleavage  planes  will  provide  missing  information  on  the  mechanism  of  slip 
and  the  low  fracture  strength.  The  atomistic  information  will  be  described  in  a  continuum 
context  developed  by  Rice  (2)  to  evaluate  the  competition  between  cleavage  fracture  and 
dislocation  nucleation  at  crack  tips.  The  model  uses  an  interplanar  potential  to  relate  the 
continuum  shear  stress  to  the  atomic  displacement  as  the  boundary  condition  to  be  satisfied  in  an 
elastic  solution  of  a  traction  free  crack.  Rice’s  analysis  leads  to  a  new  material  parameter,  the 
unstable  stacking  fault  energy,  which  is  the  maximum  energy  encountered  for  sliding  blocks  of 
material  past  one  another  on  a  slip  plane.  The  unstable  stacking  fault  energy  cannot  be  directly 
measured  by  experiment,  but  is  derivable  from  atomistic  simulations  of  stacking  faults.  Rice’s 
model  will  provide  a  mechanistic  description  for  RDX’s  observed  brittle  nature. 

The  anisotropic  material  properties  and  unstable  stacking  fault  energies  for  RDX  will  be 
determined  from  molecular  dynamics  simulations  using  a  potential  energy  function  developed  by 
Smith  and  Bharadwaj  (SB)  (6)  for  similar  nitramines  such  as  FIMX.  The  SB  potential  describes 
intramolecular  bonded  atomic  interactions  and  nonbonded  interactions  for  repulsion/dispersion 
and  electrostatics.  The  potential  is  validated  by  reproducing  experimental  RDX  material 
properties  and  structural  information  for  the  ambient  condition  a-RDX  phase  and  the  high 
pressure  y-RDX  phase.  Aided  by  high  performance  computing,  these  calculations  will  yield  for 
the  first  time  categorical  information  about  the  slip  systems  likely  to  be  activated  during  plastic 
deformation,  thereby  providing  basic  material  understanding  of  processes  that  may  lead  to  hot 
spot  formation  and  initiation  of  RDX. 


2.  Calculations 


Smith  et  al.  ( 6)  used  quantum  mechanics  calculations  to  parameterize  nitramide  and 
dimethylnitramine  (DMNA)  shown  in  figure  1(a)  and  (b)  to  develop  classical  intramolecular 
bond  potentials  commonly  used  in  molecular  dynamics.  Nitramide  and  DMNA  are  similar  in 
structure  to  the  nitroamino  functional  group  shown  in  figure  1(a)  that  makes  up  both  FIMX  and 
RDX.  Due  to  these  similarities,  a  flexible  potential  developed  for  DMNA  provides  a  building 
block  to  a  flexible  molecular  potential  for  the  more  complex  RDX  and  HMX  molecules. 
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The  quantum  mechanics  data  for  DMNA  was  used  to  fit  a  bonded  potential  energy  function 
containing  harmonic  bond  forces,  harmonic  bending  forces,  cosine  series  torsion/dihedral  forces, 
and  out-of-plane  bending  forces  of  the  O-N-O  group,  with  the  N-N  bond  modeled  as  a  harmonic 
improper  dihedral.  The  improper  dihedral  relates  the  motion  of  the  two  nitrogens  by  harmonic 
torsion  about  an  axis  between  the  two  oxygen  atoms,  which  works  to  keep  the  N-N  bond  in¬ 
plane  with  the  O-N-O  group.  The  nonbonded  atomic  interactions  were  modeled  with  an 
exponential-6  repulsion/dispersion  model,  with  parameters  given  in  other  work.  The  partial 
charges  for  Coulombic  interactions  were  found  by  fitting  an  electrostatic  potential  at  a  grid  of 
point  values  from  the  quantum  mechanics  simulations.  Simulations  using  the  DMNA  potential 
are  shown  to  reproduce  experimental  gas  phase  peaks  on  the  radial  distribution  function  for 
atomic  pairs  and  liquid  phase  vapor  pressure,  volume,  and  temperature  properties.  The  liquid 
phase  simulations  required  the  partial  charges  to  be  increased  by  25%  to  reproduce  experimental 
results. 

SB  ( 1 )  assumed  that  all  valence  and  dispersion/repulsion  terms  found  for  the  much  simpler 
DMNA  molecule  were  directly  transferable  to  the  potential  used  for  HMX.  Additional  valence 
parameters  were  then  determined  for  HMX  bonds  that  do  not  occur  in  DMNA,  namely  the  N-C- 
N  bend  and  C-N-C-N  dihedral  angles  of  the  ring.  The  potential  includes  nonbonded  coulombic 
and  dispersion/repulsion  interactions  between  all  atoms  separated  by  three  or  more  bonds, 
including  those  connected  by  1-4  dihedrals.  Partial  charges  for  the  HMX  model  were  refit  to 
reproduce  a  grid  of  quantum  mechanics  data  points  for  HMX,  with  the  constraint  that  “like” 
atoms  have  equal  charges.  The  SB  flexible  molecule  potential  was  used  to  find  the 
experimentally  unavailable  thermal  conductivity,  shear  viscosity,  and  self-diffusion  coefficients 
for  the  liquid  HMX  temperature  domain  550  <  T  <800  K,  which  are  all  important  material 
properties  for  large  length  and  time  scale  constitutive  models  used  in  continuum  hydrodynamic 
codes  (7).  The  SB  potential  was  then  used  to  model  crystalline  HMX  in  the  p-,  a-,  and  8-phases 
stable  at  ambient  conditions,  which  differ  by  molecular  packing  and  molecule  conformation  (5). 
Bedrov  et  al.  give  the  full  potential,  U,  used  for  crystalline  HMX  as 

U  =  0.5  Kfj(nj  -  rtf  +  O.SKtfdi]k  -  dtf2  +  0.5ifJfci[l  -  cos(ncpijkl )] 

+0-5 KPkl(pl-kl  +  Aijexp^-Bijrij)  -^  +  ketf  (1) 

ij  i] 
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with  parameters  listed  in  table  1 .  This  exact  form  of  the  SB  potential  and  parameter  set  from 
table  1  is  used  to  model  crystalline  RDX  in  this  work. 

Table  1.  Smith  and  Bharadwaj  potential  parameters  for  HMX/RDX. 


Bond  stretches,  U  = 

CN 

1 

LO 

O 

Bond  type 

Kfj  (kcal/mol/A2) 

rij  (A) 

O-N 

1990.1 

1.23 

N-N 

991.7 

1.36 

N-C 

672.1 

1.44 

C-H 

641.6 

1.09 

Valence  Bends,  U  =  i 

o.5/^k(et/k  -  0(°/fe)2 

Bend  type 

K®k  (kcal/mol/rad2) 

Oijk  (rad) 

O-N-O 

125.0 

2.1104 

O-N-N 

125.0 

1.8754 

N-N-C 

130.0 

1.6723 

C-N-C 

70.0 

1.843 

N-C-H 

86.4 

1.8676 

H-C-H 

77.0 

1.8938 

N-C-N 

70.0 

1.9289 

Torsions,  U  =  0.5kffei[l  —  cos  (n0i;fei)] 

Torsion  Type 

Kmi  (kcal/mol) 

n 

O-N-N-C 

8.45 

2 

O-N-N-C 

0.79 

4 

O-N-N-C 

0.004 

8 

H-C-N-C 

-0.16 

3 

C-N-C-N 

3.30 

1 

C-N-C-N 

-1.61 

2 

C-N-C-N 

0.11 

3 

Out  of  plane  bends  U  =  0.5K[jk(pfikl 

Out-of-plane  bend 
type 

K[jk  (kcal/mol/rad2) 

C-N-C...  *N 

O-N-O...  *N 

8.0 

89.3 

Where  . .  ,*N  is  the  atom  kept  in-plane 

van- 

-der-Waals  interactions,  U 

=  Al}  exp(-5t/ry)  - 

Cii/rf 

Atoms  pair  type 

Ay  (kcal/mol) 

Cy  (kcal/mol  A6) 

C-C 

14976.0 

3.090 

640.8 

C-H 

4320.0 

3.415 

138.2 

C-N 

30183.57 

3.435 

566.03 

c-o 

33702.4 

3.576 

505.6 

H-H 

2649.7 

3.740 

27.4 

H-N 

12695.88 

3.760 

116.96 

H-O 

14175.97 

3.901 

104.46 

N-N 

60833.9 

3.780 

500.0 

14 


N-0 

0-0 


67925.95 

75844.8 


3.921 

4.063 


446.6 

398.9 


Atomic  partial  charges 


Atom  type 

0 

C 

-0.5400 

N(amine) 

0.056375 

N(nitro) 

0.860625 

0 

-.458500 

H 

0.270000 

Cawkwell  et  al.  (9)  demonstrate  the  transferability  of  the  SB  potential  given  in  table  1  by 
applying  it  directly  to  simulations  of  a-RDX  under  shock  loading.  The  ambient  condition 
thermomechanical  properties,  like  the  coefficient  of  thermal  expansion,  bulk  modulus,  and 
elastic  constants  for  a-RDX,  are  not  calculated,  but  the  ambient  lattice  constants  are  reported  and 
shown  to  be  within  1.7%  of  experimental.  The  simulations  use  a  shock  front  absorbing  boundary 
layer  to  model  viscous  heating  in  shear  bands,  molecule  structure,  and  crystal  order  in  the 
amorphous  shear  band  region.  The  simulations  of  Cawkwell  et  al.  are  compared  to  RDX 
simulations  by  Bedrov  et  al.  (10)  using  the  SB  potential  with  a  modified  thermostat  and  barostat 
that  enforce  the  Hugoniot  shock  conditions  of  a  stationary  shock  wave  for  an  applied  uniaxial 
pressure. 

The  molecular  dynamics  package  DL-Poly  2.20  (11)  was  used  in  this  work  to  test  the  application 
of  the  SB  Flexible  Molecular  potential  (1)  for  the  crystalline  phases  of  a-RDX.  We  performed 
simulations  to  check  that  the  potential  recreated  the  proper  crystal  structure  of  a-RDX  by 
comparing  time  averaged  atomic  position  data  from  the  simulations  to  experimental  structural 
data  given  by  Choi  and  Prince  (12).  The  initial  atomic  coordinates  for  the  a-RDX  simulations 
were  found  by  applying  the  Pbca  space  group  symmetry  operations  to  the  asymmetric  molecule 
coordinates  given  by  Choi  and  Prince.  The  initial  configurations  for  the  y-RDX  simulations  are 
those  given  by  Davidson  et  al.  for  y-RDX  unit  cell  (13).  The  unit  cell  was  then  copied  and 
translated  to  create  a  2x3x3  crystal  lattice  in  the  (a,  b,  c)  directions,  respectively.  The  2x3x3 
RDX  crystal  contains  18  unit  cells,  144  molecules,  and  3,024  atoms.  Parallelepiped  periodic 
boundary  conditions  are  used  in  the  simulations,  which  allow  for  non-orthogonal  lattice  vectors 
of  different  lengths.  Periodic  boundary  conditions  are  used  to  simulate  an  infinite  crystal  in  all 
directions  in  order  to  get  bulk  thermodynamic  properties  free  from  the  effects  of  surfaces. 

The  real  space  cut-off  for  the  non-bonded  Van  der  Waals  and  electrostatic  interactions  was  set  to 
10  A,  which  is  large  enough  to  allow  for  interactions  between  nearest  neighbor  molecules.  The 
smallest  simulation  cell  dimension  must  be  at  least  twice  as  large  as  the  10  A  real  space  cut-off 
requiring  a  2x3x3  unit  cell  crystal.  In  DL-POLY,  a  long  range  correction  to  the  potential  energy 
is  applied  to  the  short  range  van  der  Waals  interactions  to  account  for  the  attractive  forces 
between  atoms  at  distances  greater  than  the  real  space  cut-off.  The  long-range  electrostatic 
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interactions  are  calculated  using  the  Ewald  sum  method,  which  splits  the  potential  up  into  two 
parts:  atoms  separated  by  a  distance  less  than  the  cut-off  are  treated  using  a  direct  Coulomb  sum, 
and  atoms  separated  at  distances  larger  than  the  cut-off  are  treated  using  a  Fourier  series  sum  in 
reciprocal  space.  In  this  work,  DL-POLY  is  used  to  automatically  set  the  Ewald  parameters  by 
specifying  the  precision  of  the  relative  error  in  the  convergence  of  the  real  space  sum  to  0.3e-6. 
DL-POLY  uses  the  relative  error  to  automatically  optimize  the  reciprocal  space  k  vectors  and  the 
Ewald  convergence  parameter.  In  these  simulations,  the  Ewald  sum  is  shown  to  be  adequately 
converged  by  checking  that  the  Coulombic  virial  is  the  negative  of  the  Coulombic  energy. 

The  simulations  were  run  using  the  verlet  leapfrog  time  integrator.  The  highest  vibrational 
frequency  in  the  crystal  is  the  C-H  bond  stretch,  which  has  a  period  of  approximately  12fs.  In 
order  to  conserve  the  total  system  energy,  approximately  1 0  integration  steps  should  be  taken  per 
period,  and  in  this  work,  a  0.75-fs  timestep  was  used.  This  integration  step  could  have  been 
increased  by  treating  the  C-H  bond  as  rigid,  which  may  be  done  in  future  work. 

All  simulations  in  this  work  use  the  Nose-Hoover  thermostat  to  control  the  simulation 
temperature.  The  Nose-Hoover  thermostat  functions  by  scaling  the  velocity  using  a  scalar 
friction- like  term  that  is  coupled  to  an  outside  heat  bath.  For  these  simulations  the  temperature  is 
set  to  ZK500K.  The  coupling  time  constant  for  the  thermostat  in  these  simulations  is  1.0  ps, 
which  is  related  to  the  thermostat’s  “effective”  mass. 

This  work  contains  both  constant  pressure/temperature  (NPT)  and  constant  volume/temperature 
(NVT)  simulations.  In  the  constant  pressure  simulations,  the  pressure  is  controlled  with  a 
barostat  coupled  to  the  thermostat  through  a  modification  of  the  equations  of  motion.  The 
barostat  controls  the  size  and  shape  of  the  simulation  cell  to  maintain  the  prescribed  pressure. 

The  barostat  is  applied  as  a  friction-like  term  that  scales  the  velocity  and  simulation  cell  size  and 
shape.  In  an  isobaric-isothermal  simulation  (NPT),  the  pressure  is  controlled  by  a  scalar  friction 
term  scaling  the  cell  shape.  This  work  uses  an  iso-stress  barostat  (NST),  which  uses  a  second 
order  tensor  friction  term  that  scales  the  shape  and  size  of  the  simulation  cell  independently.  The 
NST  with  the  parallelepiped  periodic  boundary  conditions  allows  for  the  simulation  to  proceed 
under  the  least  geometrically  constrained  boundary  conditions,  and  to  reach  an  equilibrated 
structure  predicted  by  the  molecular  potential  at  the  applied  temperature  and  pressure.  The 
barostat  coupling  time  in  these  simulations  is  set  to  1.0  ps,  which  is  related  to  the  barostat’s 
“effective”  mass. 

All  simulations  are  started  with  an  initial  7.5  ps  equilibration  period,  where  the  random  velocities 
used  to  seed  the  simulation  are  scaled  every  five  integration  steps  to  maintain  the  correct 
temperature.  Total  system  energy  is  not  conserved  during  this  part  of  the  simulation.  After 
equilibration,  the  temperature  scaling  is  turned  off  and  the  system  evolves  according  to  the 
equations  of  motions  in  their  respective  ensembles.  This  period  is  referred  to  as  thermalizing  in 
this  work,  and  during  this  period  the  ensemble  energy  is  conserved,  and  meaningful  time 
averaged  values  and  standard  deviations  are  collected. 
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We  performed  two  separate  sets  of  NST  simulations,  starting  at  two  different  states  to  study 
effect  of  pressure  on  RDX.  The  first  set  started  in  the  ambient  a-RDX  crystal  configuration 
given  by  Choi  and  Prince  (12)  at  P=1  atm  and  T=300  K,  and  the  pressure  is  incrementally 
increased  from  0.0-5 .2  GPa  in  increments  of  0.5GPa.  The  second  set  of  simulations  starts  in  the 
y-RDX  configuration  given  by  Davidson  et  al.  (13)  at  P=5.2  GPa  and  T=300  K,  and  the  pressure 
is  incrementally  decreased  from  5. 2-0.0  GPa  in  increments  of  0.5  GPa.  The  final  atomic 
configuration  and  lattice  vectors  from  the  subsequent  pressure  simulation  are  used  as  the  initial 
conditions  at  the  new  pressure.  Each  simulation  is  equilibrated  for  7.5  ps  and  thermalized  for 
100  ps. 

Isostrain  NVT  simulations  were  also  performed  here  to  study  the  orientational  dependence  of  the 
system  state  to  applied  uniaxial  strain.  The  system  was  strained  by  applying  displacements  to 
molecule  centers  of  mass  and  lattice  vectors  according  to  the  Lagrange  strain  description.  The 
strain  components,  En,  E22  and  E33,  were  incremented  individually,  requiring  three  separate  sets 
of  simulations.  The  thermalized  average  lattice  vectors  and  atomic  configuration  from  the  NST 
simulations  for  T=300  K  and  P=0  for  a-RDX,  or  P=5.2  GPa  for  y-RDX,  were  used  as  the  initial 
zero-strained  configuration.  The  strain  was  incrementally  increased  from  the  subsequent  steps 
final  atomic  configuration  in  strain  increments  of  1%  up  to  10%  strain.  Each  simulation  is 
equilibrated  for  7.5  ps  and  thermalized  for  75  ps. 

A  third  set  of  NVT  annealing  simulations  were  done  to  study  the  effect  of  atomic  disregistry  or 
stacking  mismatch  on  a-RDX  slip  planes.  During  the  NVT  annealing  process,  the  system  was 
thermalized  to  25  K  for  100  timesteps,  after  which  the  force  is  minimized  using  conjugate 
gradient,  and  then  re-thermalized  and  minimized.  This  thermalize-minimize  process  is  done  20 
times,  and  the  minimum  energy  configuration  of  these  cycles  is  used  for  postprocessing.  The 
initial  configurations  for  these  simulations  come  from  annealing  a  3x3x3  a-RDX  system  with 
initial  coordinates  from  Choi  and  Prince  (12).  The  3x3x3  a-RDX  lattice  is  then  used  to  create  a 
sandwich  structure  containing  30  unit  cells  in  the  direction  normal  to  the  slip  plane,  with  the 
atomic  positions  of  the  top  14  and  bottom  14  unit  cells  frozen  in  place.  A  4  mn  vaccum  layer  is 
added  to  the  top  of  the  30  unit  cells.  The  atomic  positions  of  the  unit  cells  above  the  slip  plane 
are  then  shifted  by  increments  of  1/1 0th  of  a  unit  cell  to  create  a  grid  of  100  configurations 
shifted  in  the  slip  plane.  The  simulation  uses  three-dimensional  (3-D)  orthorhombic  periodic 
boundary  conditions.  The  single  layers  of  flexible  cells  above  and  below  the  slip  plane  are 
allowed  to  rotate  and  translate,  and  their  intramolecular  degrees  of  freedom  are  allowed  to  relax 
during  the  thermalization  and  minimization  procedure,  hopefully  not  getting  trapped  in  a  local 
minimum  and  relaxing  to  the  global  minimum.  The  frozen  cells  are  used  simulate  bulk  material 
by  shielding  the  flexible  unit  cells  from  the  free  surfaces,  and  also  to  constrain  the  molecule 
center  of  mass  motion  to  be  mostly  normal  to  the  slip  plane. 
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3.  Results  and  Discussion 


3.1  NST  Pressurization  Simulations 


The  NST  simulations  in  this  work  are  used  to  study  the  effect  of  pressure  on  the  bulk  crystal  at 
constant  temperature.  Experiments  and  simulations  run  under  these  conditions  are  used  to  find 
the  bulk  modulus,  B,  of  the  material  given  by 


B  =  -V 


(2) 


where  V  is  the  volume,  P  is  the  pressure,  and  T  is  the  temperature.  The  bulk  modulus  is  an 
intensive  thermodynamic  material  property  that  provides  a  scalar  measure  of  the  crystal’s 
volumetric  response  to  pressure.  In  this  work,  the  isothermal  pressure  versus  volume  data  from 
simulation,  shown  in  figure  2,  is  fit  to  the  third  order  Birch-Mumaghan  equation  of  state 
(BMEOS)  to  find  the  bulk  modulus  and  its  derivative  (14)  given  by 

7/  3  /14A5/3l  (\  .  3  ^\fVn\~2/  3 


P  00  = 


3  B0 


©  "  -  © 


1  +  “  (.B'0  —  4) 


-  1 


(3) 


where  V0  is  the  volume  at  zero  pressure,  and  B„  and  B'0  are  the  bulk  modulus  and  its  derivative  at 
zero  pressure.  This  method  gives  the  bulk  modulus  at  zero  pressure  and  is  only  applied  to  a- 
RDX  data  shown  by  the  filled  symbols.  Values  for  B0  and  B'0  are  given  in  table  2,  and  are  shown 
to  be  in  agreement  with  experimental  values. 


Figure  2.  Pressure  versus  normalized  volume.  Circles  are  simulation  results,  triangles  are  experimental 
results  experiement  (75).  Filled  symbols  are  for  a-RDX  and  unfilled  symbols  are  for  y-RDX. 
Lines  connecting  simulation  data  points  are  only  to  highlight  the  phase  transition  for  y-  to 
a-RDX. 
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Table  2.  Birch-Murnaghan  equation  of  state  parameters. 


q-RDX  3rd  order  BMEOS  parameters 


P  range  (Gpa) 

B0(GPa) 

B0' 

V0(A3) 

This  work 

0-15 

13.0 

9.2 

1634.4 

Olinger3 

0-3.95 

12.1 

8.6 

1641 

Davidson*3 

9.8 

11.4 

1614.1 

A  phase  transition  from  a-  to  y-RDX  is  reported  by  Davidson  et  al.  (13)  to  occur  at  3.9  GPa. 

The  experimental  isothermal  compression  data  of  Olinger  et  al.  (15),  highlighted  by  the  triangle 
data  points  in  figure  2,  show  the  phase  transition  to  occur  between  4  and  5  GPa.  Two  separate 
curves  develop  for  the  simulation  data  in  figure  2  for  pressures  above  2GPa.  The  compression 
simulation  of  a-RDX,  shown  by  the  solid  line  and  filled  circles,  does  not  show  a  kink  in  the 
volumetric  compression  data,  and  a  phase  transition  is  likely  not  occuring.  A  small  kink 
develops  in  the  decompression  simulation  of  y-RDX,  shown  by  the  converging  of  the  hollow 
circles  and  dashed  line  with  the  solid  line  of  a-RDX,  between  2.5  and  2  GPa.  The  kink  signifies 
the  phase  transition  between  the  thennalized  y-  and  a-RDX  phases  found  from  simulation.  A 
much  larger  volume  change  occurs  experimentally  than  that  observed  from  the  simulations. 

The  accuracy  with  which  the  SB  potential  predicts  the  lattice  constants  for  the  a-  and  y-RDX 
phases  is  presented  in  figure  3  and  tabulated  in  table  3.  The  point  of  comparison  for  a-RDX  is 
the  ambient  state  at  P=0  and  T=300  K,  which  was  the  state  measured  experimentally  by  Choi  and 
Prince  (12).  The  y-RDX  point  of  comparison  is  P=5.2  and  T=300  K,  the  experimental  conditions 
for  the  reported  y-RDX  structure  of  Davidson  et  al.  (13). 


Figure  3.  Lattice  constants  versus  pressure.  Diamonds  (a),  circles  (b),  and  squares  (c)  are  results  from 

simulations;  triangles  are  from  experiment  (75).  Filled  symbols  are  for  a-RDX  and  unfilled  symbols 
are  for  y-RDX.  Lines  connecting  simulation  data  points  are  only  to  highlight  the  phase  transition  for 
y-  to  a-RDX. 
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Table  3.  Average  lattice  constants  and  volume.  Percent  error  compared  to  experiment  in  parentheses. 


Lattice  Dimensions 


<a>  (A) 

<b>  (A) 

<c>  (A) 

<V>  (A3) 

-Exptr 

13.182 

11.574 

10.709 

1633.8 

-This  work 

13.446  (2.00%) 

11.531  (-0.37%) 

10.534  (-1.63%) 

1633.3  (-0.03%) 

SAPTb(298) 

13.259  (0.58%) 

11.634  (0.52%) 

10.754  (0.42%) 

1658.9  (1.54%) 

-Exptf 

12.565 

10.930 

9.477 

1301.5 

-This  work 

12.685  (1.0%) 

11.063  (1.2%) 

9.644  (-1.8%) 

1353.3  (4.0%) 

aChoi,  Prince  (72)  bPodesczwa,  Rice,  Szalewicz  (16)  cDavidson  et  al.  (13) 


The  a-RDX  lattice  constants,  shown  by  the  filled  symbols  and  solid  line,  are  overpredicted  in  the 
a-direction,  but,  volumetrically,  this  is  balanced  out  by  the  underprediction  in  the  c-direction. 
There  is  also  not  a  sharp  kink  in  the  a-RDX  lattice  constants,  signifying  the  lack  of  a  phase 
transition.  The  phase  transition  in  y-RDX  simulations,  shown  by  the  unfilled  symbols  and 
dashed  lines,  is  apparent  in  the  lattice  constants  by  the  divergence  of  the  dashed  and  solid  lines  at 
P  >  2  GPa.  The  y-RDX  phase  more  closely  matches  the  experimental  data  points  than  the  a- 
RDX  simulations.  Across  the  phase  transition,  the  y-RDX  simulations  show  a  decrease  in  the  c- 
direction,  an  increase  in  the  b-direction,  and  little  change  in  the  a-direction,  all  of  which  are 
observed  experimentally.  The  inability  of  the  SB  potential  to  capture  the  large  volume  change 
across  the  phase  transition  in  figure  2  is  caused  by  the  errors  in  a-RDX  lattice  constants. 

The  Radial  Distribution  Function  (RDF)  is  also  used  here  to  detennine  the  phase  of  the 
simulations  due  to  uncertainty  seen  in  volumetric  data  in  figure  2  and  lattice  constants  in  figure  3 
for  the  a-  and  y-RDX  phases.  The  RDF  is  used  to  plot  the  molecular  density  as  a  function  of 
distance  from  a  particular  atom.  The  RDF  is  scaled  by  the  molecule  density  and,  so,  will 
approach  unity  at  large  distances.  Spikes  in  the  RDF  signify  atoms.  The  RDF  in  this  work  was 
created  from  30  realizations  of  the  configuration  during  the  thermalization.  The  experimental 
RDFs  shown  by  the  dashed  lines  in  figure  4  show  sharp  peak  because  the  molecules  lack  thermal 
vibrations.  Thermalizations  of  the  states  that  the  a-  and  y-RDX  crystal  structure  were 
experimentally  determined  at  are  shown  by  the  solid  colored  line,  T=300K,  P=0  for  a-RDX;  by 
the  solid  blue  line  at  the  bottom;  and,  P=5.2GPa  for  y-RDX,  shown  by  the  solid  red  line  at  the 
top.  Figure  4  shows  the  RDFs  for  the  pressure  data  points  plotted  in  figures  2  and  3  for  the 
pressurization  of  a-RDX.  The  pressure  of  each  RDF  is  labeled  in  figure  4. 
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Figure  4.  Scaled  molecule  center  of  mass  RDF  for  compression  simulations  of  a-RDX.  Each  RDF  is  at  a 
different  pressure  shown  by  labels.  Top  and  bottom  thermalized  (solid  line)  and  experimental 
(dashed  line)  (12)  RDF’s  in  blue  (a-RDX)  and  red  (y-RDX)  shown  for  reference. 

The  experimental  RDF  of  a-RDX  is  compared  to  its  thermalized  RDF  at  the  same  state,  T=300 
K,  P=1  atm,  shown  by  the  blue  dashed  and  solid  line,  respectively.  The  first  four  peaks  of  the 
experimental  data  are  accurately  captured  by  the  thermalized  data.  Two  distinct  peaks  in  the  8- 
10  A  region  join  into  one  broad  peak  in  the  simulation  data.  Molecules  within  10  A  make  up  the 
first  set  of  nearest  neighbors,  and  this  distance  is  specified  as  the  real  space  cut-off  of  nonbonded 
interactions  in  the  simulation  set-up.  Beyond  about  12  A,  the  RDF  approaches  unity.  The 
sequence  of  pressurization  simulations  for  a-RDX,  shown  in  figure  4,  start  in  the  blue 
configuration  of  a-RDX,  shown  at  the  bottom  of  the  plot.  The  solid  blue  line  is  the  same  RDF  as 
the  black  line  labeled  P=0.0.  The  pressure  is  then  increased,  and  the  peaks  of  the  RDFs  move 
closer  to  the  origin  as  the  molecules  are  squeezed  closer  together.  During  pressurization,  peaks 
narrow  and  separate  as  the  thermal  vibrations  are  reduced.  At  around  P=4  GPa,  the  peaks  begin 
to  sharpen  in  the  5-8  A  region,  and  the  broad  peak  from  8-10  A  narrows,  but  this  is  result  of  the 
reduced  molecular  mobility  as  the  pressure  is  increased,  and  not  a  phase  transition.  The  top  solid 
black  line  labeled  P=5.2  GPa  is  not  representative  of  the  of  the  solid  red  line  of  the  thermalized 
y-RDX  phase  at  5.2  GPa.  There  is  a  difference  in  the  number  of  peaks  in  the  8-10  A  region 
(3  for  y-RDX,  1  for  a-RDX),  the  relative  size  of  peaks  from  6-10  A  (first  peak  is  largest  for  y- 
RDX,  first  peak  is  smallest  for  a-RDX),  and  the  location  of  the  first  peak  at  4  A. 

The  RDF  for  g-RDX  decompression  with  the  same  reference  state  for  a-  and  g-RDX  in  blue  and 
red,  respectively,  is  shown  in  figure  5.  In  these  simulations  the  pressure  is  reduced  and  starts  in 
the  configuration  labeled  P=5.2.  As  the  pressure  is  reduced,  a  noticeable  change  occurs  at  P=2.0 
GPa,  which  is  also  the  pressure  that  kinks  in  the  volume,  and  lattice  constants  were  seen  in 
figures  2  and  3.  For  P  <  2  GPa,  the  system  goes  to  the  thermalized  a-RDX  configuration. 
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Figure  5.  Scaled  molecular  center  of  mass  RDF  for  decompression  simulations  of  y-RDX.  Each  RDF  is  at  a 
different  pressure  shown  by  labels.  Top  and  bottom  thermalized  (solid  line)  and  experimental 
(dashed  line)  (75)  RDF’s  in  blue  (a-RDX)  and  red  (y-RDX)  shown  for  reference. 

A  phase  transition  determined  by  the  locations  of  peaks  in  the  RDF  can  be  summarized  by  a 
single  data  point  determined  by  the  time-averaged  mean  square  displacement  (MSD)  of  atoms. 
The  MSD  compares  the  current  atomic  location  to  the  initial  atomic  location  at  the  start  of  the 
simulation.  In  this  work,  the  MSD  is  averaged  for  each  atomic  species — carbon  (C),  nitrogen 
(N),  oxygen  (O),  and  hydrogen  (FI).  The  MSD  of  the  carbon  atoms  contained  in  the  ring  is 
plotted  in  figure  6  and  provides  a  good  approximation  to  the  MSD  of  the  molecule  center  of 
mass.  Since  each  set  of  compression  or  decompression  simulations  is  done  incrementally,  a 
large  spike  in  the  MSD  will  occur  when  the  molecules  shift  from  the  phase  they  are  started  in  to 
a  new  one,  as  is  shown  for  the  y-RDX  simulation  data  at  P=2  GPa.  The  time-averaged  MSD  for 
a-RDX  does  not  show  a  spike  in  figure  6,  and  it  can  be  concluded  that  a  phase  transition  does 
not  occur  during  the  pressurization  of  a-RDX  in  these  simulations,  even  though  the  reverse 
transition  from  y-  to  a-RDX  readily  occurs  at  ~2  GPa.  The  difference  between  the  a-  and  y- 
RDX  MSD  at  pressures  below  2  GPa  could  be  caused  by  a  rotation  in  the  molecule.  The  MSD  is 
a  measure  of  the  fluctuations  in  atomic  positions.  Other  fluctuation  formulas  exist  for  energy 
terms  that  can  be  used  to  determine  thermodynamic  properties,  such  as  the  heat  capacity,  which 
can  also  be  used  to  provide  evidence  of  a  phase  transition. 
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Pressure  (Gpa) 


Figure  6.  Time  averaged  MSD  for  a-  and  g-RDX  versus  pressure. 


Up  to  this  point,  the  determination  of  the  crystal  structure  has  been  based  on  intermolecular 
properties,  such  as  molecule-to-molecule  spacing,  lattice  constants,  and  volume.  None  of  these 
measures  highlight  the  main  difference  between  the  a-  and  y-RDX  crystals,  which  is  contained 
in  the  molecule  conformation  and  the  orientation  of  the  nitro  groups  with  respect  to  the  central 
amine  ring.  This  orientation  is  referred  to  as  the  wag  angle  and  is  shown  below  in  figure  7.  The 
orientation  of  the  N-N  bond  with  the  C-N-C  plane  is  the  wag  angle,  5,  shown  in  figure  7.  Wag 
angles  in  this  work  are  referenced  by  their  respective  nitro  groups  based  on  the  numbering  shown 
for  a  right  handed  molecule  in  figure  7(b). 


Figure  7.  (a)  Wag  angle,  5,  used  to  describe  nitro  group  orientation,  (b)  N4-N7  wag  angle  shown  on 
thermalized  right  handed  molecule. 

RDX  molecules  are  described  by  8  for  each  N-N  group — Equatorial  (E)  for  8<0°,  Axial  (A)  for 
8>10°,  and  Intermediate  (I)  for  0°<8<10°.  These  definitions  of  A,  E,  and  I  are  based  on 
somewhat  arbitrary  values  of  8  that  only  apply  to  RDX.  Choi  and  Prince  (72)  determined  the  a- 
RDX  unit  cell  to  belong  to  the  Pbca  space  group  and  contain  eight  symmetrically  equivalent 
molecules  with  nitro  groups  in  the  Axial,  Axial,  Equatorial  (AAE)  positions.  Davidson  et  al. 
(73)  determined  the  y-RDX  unit  cell  to  belong  to  the  Pca2i  space  group,  with  four  symmetrically 
equivalent  sites,  and  each  site  containing  two  separate  molecules  in  two  different  configurations, 
AAE  and  AAI. 
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In  this  work,  instantaneous  wag  angles  are  determined  for  every  N-N  group  in  the  simulation  cell 
at  30  different  instantaneous  atomic  configurations  separated  by  0.75ps.  Each  wag  angle  type  is 
then  summed  according  to  which  of  the  eight  molecule  types  of  the  unit  cell  it  belongs  to. 
Histogram  data  of  each  wag  angle  type  for  each  of  the  eight  molecules  of  the  unit  cell  is  then 
plotted  to  determine  the  conformation  of  the  molecule,  as  is  shown  in  figures  8(a)-(c)  for 
a-RDX,  and  in  figures  9(a)-(c)  for  y-RDX.  Only  the  outline  of  the  histogram  data  is  plotted  in 
figure  8. 


Figure  8.  a-RDX  AAE  conformation  determined  from  Wag  angle  histograms  at  T=300K,  P=0,  for  (a)  N4-N7, 
(b)  N5-N8,  and  (c)  N6-N9  for  each  molecule  type  in  the  unit  cell. 


Figure  9.  y-RDX  AAE  &  AAI  conformations  determined  from  Wag  angle  histograms  at  T=300K,  P=5.2GPa, 
for  (a)  N4-N7,  (b)  N5-N8,  and  (c)  N6-N9  for  each  molecule  type  in  the  unit  cell. 

The  difference  between  the  a-  and  y-RDX  phase  is  highlighted  by  the  separation  of  the 
histogram  data  for  the  N4-N7  and  N6-N9  wag  angles  in  figure  9(a)  and  (c),  signifying  two 
distinct  molecule  conformations  in  the  unit  cell.  The  peak  thermalized/experimental  wag  angles 
for  a-RDX  are  25°/34°,  45°/33°,  and  -20°/-20°  for  the  N5-N8  axial  (A),  N6-N9  axial  (A),  and 
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N4-N7  equatorial  (E),  respectively.  The  peak  thermalized/experimental  wag  angles  for  type  1 
molecules  of  y-RDX  are  36°/36°,  40°/36°,  and  8°/10°  for  the  N5-N8  axial  (A),  N6-N9  axial  (A), 
and  N4-N7  equatorial  (I),  respectively.  The  peak  thermalized/experimental  wag  angles  for  type 
2  molecules  of  y-RDX  are  38°/40°,  20°/17°,  and  —1 4°/-2°  for  the  N5-N8  axial  (A),  N6-N9 
axial  (A),  and  N4-N7  equatorial  (E),  respectively.  The  unit  cell  of  for  a-RDX  is  shown  in 
figure  10(a)-(c),  and  using  the  same  basis  molecules,  the  unit  cell  of  y-RDX  is  shown  in 
figure  10(d)—  (f). 


Figure  10.  RDX  unit  cells  with  H  and  O  atoms  removed  projected  onto  crystal  planes  (shown  in  parentheses). 
Molecules  colored  according  to  y-RDX  molecule  type.  Type  1  AAI  molecules  shown  in  blue  and 
type  2  AAE  molecules  shown  in  red. 

The  AAE  molecules  shown  in  red  undergo  very  little  change  in  lattice  location  or  confonnation, 
as  shown  by  comparing  the  a-  to  y-RDX  unit  cells  in  figure  10.  The  change  between  an  AAE 
and  AAI  molecule  is  shown  by  comparing  the  blue  molecules  in  figures  10(a)  and  (d).  The 
rotation  required  to  allow  for  the  AAI  confonnation  is  shown  by  comparing  the  blue  molecules 
in  figures  10(b)  and  (e).  The  AAI  or  type  2  molecules  undergo  a  -20°  rotation  about  an  axis 
oriented  normal  to  the  amine  ring  when  compared  to  the  orientation  in  a-RDX.  This  molecule 
rotation  seems  to  be  essential  to  the  phase  transition,  and  isotropically  stressing  the  a-RDX 
simulation  cell  may  not  allow  this  rotation  to  occur.  The  pressurization  simulations  starting  in  a- 
RDX  may  not  go  through  a  phase  transition  due  to  some  type  of  artificial  constraint  on  molecule 
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rotation  caused  by  the  simulation  cell  or  boundary  conditions.  The  a-  to  y-RDX  phase  transition 
may  be  possible  to  recreate  by  applying  non-hydrostatic  states  of  stress  or  strain,  which  could 
free  up  rotational  degrees  of  freedom  in  the  molecule. 


3.2  NVT  Uniaxial  Strain  Simulations 


The  orientation  dependence  of  the  RDX  crystal  response  to  uniaxial  loading  will  be  studied  in 
this  work  by  using  NVT  simulations  of  uniaxially  strained  RDX  lattices  in  the  a-  and  y-RDX 
phases.  We  will  determine  the  orthotropic  elastic  constants  for  linear  elasticity  for  these  phases 
and  the  dependence  of  strain  on  the  phase  transition.  The  strain  orientation  dependence  of  the 
phase  transition  will  elucidate  the  mechanism  to  which  it  occurs  and  possibly  provide  an 
explanation  of  why  it  does  not  occur  under  isotropic  NPT  pressurization.  The  uniaxial  strain 
simulations  here  will  provide  insight  into  highly  constrained  deformation  that  could  be  caused  by 
high  pressure  shock  fronts  and  rigid  boundary  conditions. 


The  initial  configurations  for  the  NVT  simulations  come  from  NST  simulations  of  the  a-  and  y- 
RDX,  thermalized  at  the  experimental  states  at  which  they  were  measured — T=300  K,  P=0  for 
a-RDX,  and  T=300  K  and  P=5.2  GPa  for  y-RDX.  The  b-  and  c-axes  for  y-RDX  have  been 
switched  from  those  given  by  Davidson  et  al.  (13)  in  order  for  them  to  correspond  to  the  a-RDX 
crystal  axes  given  by  Choi  and  Prince  (12).  The  thermalized  a-  and  y-RDX  crystals  were 
isostrained  by  each  strain  component  in  compression  and  tension,  respectively.  These  strained 
configurations  were  then  thermalized  in  the  NVT  ensemble.  Ten  strain  increments  of  0.1%  were 
applied  to  the  initial  configurations  up  to  1 .0%  stain.  The  time-averaged  stress  tensor  from  these 
simulations  was  used  to  determine  the  relationship  between  stress  and  applied  strain.  Assuming 
the  strain  increments  are  small,  and  the  material  is  orthotropic  due  to  the  symmetry  of  the  a-  and 
y-RDX  space  groups,  the  material  can  be  approximated  as  an  orthotropic  linear  elastic  material 
with  the  strain  relationship  given  in  Voight  notation  as 
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where  Qj  are  the  elastic  constants,  and  the  a-,  b-,  and  c-directions  are  given  by  the  ei,  e?,  and  e3 
Cartesian  coordinates.  Orthotropic  elastic  constants  from  these  simulations  are  given  in  table  4, 
and  compared  to  a-RDX  experimental  elastic  constants  and  other  simulation  data.  The  y-RDX 
elastic  constants  are  measured  from  P=5.2  GPa  and,  therefore,  are  much  more  compressed  and 
stiffer  than  the  a-RDX  elastic  constants.  When  atoms  are  brought  closer  together  due  to  large 
compressions,  the  atoms  repel  one  another — caused  by  Pauli  exclusion  principle — which  is 
captured  in  the  SB  potential  function  by  the  exponential  term  of  the  Buckingham  potential.  Both 
sets  of  elastic  constants  show  similar  trends,  and  the  ordering  of  stiffest  directions  stays  constant 
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between  phase,  Cn>C22>C33.  The  largest  increase  is  in  the  coupling  or  Poisson  terms  affecting 
the  extensional  strains,  C23,  C31,  and  C12.  There  is  not  experimental  data  to  verify  the  y-RDX 
elastic  constants,  but  the  a-RDX  elastic  constants  compare  well  to  other  a-RDX  simulations  by 
Sewell  et  al.  ( 1 7)  using  a  different  potential  and  the  experimental  data. 

Table  4.  Elastic  tensor:  b-  and  c-axes  of  y-RDX  have  been  switched  from 
that  given  by  Davidson  et  al.  (13)  in  order  to  correspond  to  a-RDX 
coordinates  given  by  Choi  and  Prince  (72). 


Orthotropic  elastic  constants  (GPa) 


(y)NVT 

(a)NVT 

(a)Sewell“ 

(a)Haycraftb 

(a)Schwarz0 

Cn 

80.28 

25.04 

26.87 

36.67 

25.60 

C22 

67.02 

23.80 

24.08 

25.67 

21.30 

C33 

57.85 

23.35 

17.69 

21.64 

19.00 

C44 

11.95 

3.12 

8.40 

11.99 

5.38 

C55 

16.30 

7.69 

5.30 

2.72 

4.27 

c66 

13.43 

5.18 

7.60 

7.68 

7.27 

C23 

45.56 

6.26 

6.32 

9.17 

6.40 

C31 

36.95 

7.59 

5.68 

1.67 

5.72 

C12 

37.83 

11.24 

6.27 

1.38 

8.67 

“Sewell  et  al.  (77),  bHaycraft  et  al.  (18),  cSchwarz  et  al.  (19) 


The  orthotropic  elastic  constants  are  accurately  predicted  using  uniaxially  strained  lattices  and 
the  NVT  ensemble  for  the  a-RDX  thermalized  crystal.  This  method  should  be  extendable  to 
larger  strains  where  the  stress/strain  relationship  becomes  nonlinear.  The  diagonal  strain 
components  between  a-RDX  at  P=0  and  y-RDX  at  P=5.2  GPa  are  6%,  4%,  and  10%  for  En,  E22, 
and  E33,  respectively.  The  strain  components  for  a-RDX  at  P=2GPa,  the  simulated  a-  to  y-RDX 
phase  transition  pressure,  are  -3%,  -3%,  and  -4%  for  En,  E22,  and  E33,  respectively.  From  the 
strain  found  for  the  NST  simulations,  it  is  determined  that  the  system  should  be  strained  up  to 
10%  in  compression  for  a-RDX  and  10%  extension  for  g-RDX. 

The  simulation  procedure  uses  the  final  atomic  coordinates  of  the  previous  strained  simulation 
and  deforms  them  by  1.0%.  Deformation  is  determined  using  the  deformation  gradient,  Fy, 
given  by 

g=[h][h0]-\  (5) 

where  [h]  and  [h0]  are  tensors  containing  column  wise  lattice  vectors.  The  configuration  to  be 
strained  is  given  by  [hQ]  and  the  final  configuration  for  the  new  simulation  is  [h].  The  Lagrange 
strain,  Ey,  for  this  is  given  by 

g  =  l{gTg-L)  =l([ho]-T[h]T[h][h0]-i-lJ,  (6) 
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where  I  is  the  identity  tensor.  The  system  is  incrementally  strained  so  [h0]  changes  to  the 
subsequent  strained  configuration  that  [h]  is  calculated  from.  F;j  is  taken  to  be 


1  +  a  0 
F  =  0  1  +/? 

“  .  0  0 
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1  +Y- 


(7) 


where  a,  [1,  and  y  are  applied  individually  and  are  set  to  a  constant  value  of  0. 1 .  The  new 
configuration,  [h],  is  given  by 

[h]  =  F[h0]  (8) 


The  molecule  centers  of  mass  are  displaced  according  to  and  all  atoms  of  the  molecule  are 
moved  as  a  rigid  body  so  as  not  to  strain  the  intramolecular  bonds.  Each  simulation  is 
equilibrated  for  7.5  ps  and  thermalized  for  75  ps.  Time-averaged  values  from  the  thermalization 
portion  are  used  in  this  analysis,  and  all  values  are  offset  so  that  the  values  for  a-RDX  at  P=0.0 
are  zero.  All  energies  are  given  per  mole  of  RDX  molecules.  Only  the  results  for  the  a-RDX 
uniaxial  strain  simulations  will  be  presented.  Some  difficulties  were  found  in  determining  a 
common  reference  state  between  the  a-  and  y-RDX  uniaxial  strained  states. 

Final  strained  a-RDX  unit  cells  for  E;i=-9. 1%  are  shown  in  figure  1 1 .  The  projection  axis  of 
each  figure  is  labeled  at  the  tops  of  the  columns,  and  the  strain  direction  is  labeled  to  the  left  of 
each  row.  The  top  row  is  the  unstrained,  thermalized  a-RDX  configuration.  The  second  row  is 
for  Ei  i=-9. 1%,  and  very  little  change  is  observed.  The  third  row  is  for  E22=-9. 1%,  and  the 
location  of  the  two  central  molecules  shift  in  the  b-direction  shown  for  the  (001)  projection.  The 
fourth  row  is  for  E33=-9.1%;  the  (010)  project  shows  the  AAI  and  AAE  molecule  conformations, 
and  the  (001)  project  shows  the  rotation  of  the  AAI  molecules  as  is  seen  in  y-RDX. 

For  all  of  the  cases  of  uniaxial  strain,  phase  transitions  are  observed  in  the  E22  and  E33  direction, 
and  only  the  E33  strain  produces  the  phase  transition  to  y-RDX.  Straining  in  the  E33  or  along  the 
c-direction  while  holding  the  a-  and  b-lattices  fixed  must  provide  the  required  molecule  mobility 
needed  for  the  a-  to  y-RDX  phase  transition.  Uniaxial  straining  along  the  different  crystal  axes 
activates  different  terms  of  the  SB  potential.  The  overall  value  of  the  SB  potential  energy  is 
plotted  in  figure  12  for  each  strain  direction.  The  large  grouping  of  data  points  near  the  origin  is 
for  the  linear  elastic  strain  calculations  used  to  find  the  elastic  constants  presented  in  table  4.  The 
other  large  grouping  of  points  for  the  E22  and  E33  curves  are  in  the  region  of  the  phase  transition. 
The  curvature  of  the  potential  energy  versus  strain  changes  at  the  phase  transition.  The 
magnitude  of  the  potential  energy  for  the  E22  and  E33  strain  directions  are  of  similar  value,  which 
is  not  surprising  because  they  have  similar  stiffness  values  for  C22  and  C33  in  table  4.  The 
potential  energy  for  the  En  direction  is  the  highest  and  it  is  also  the  stiffest  direction.  The  RDF 
and  wag  angles  used  to  determine  the  y-  to  a-RDX  transition  for  the  NST  decompression 
simulations  verify  that  the  phase  from  E33  strain  is  the  y-RDX  phase. 
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Figure  11.  Final  strained  a-RDX  unit  cells.  Projection  planes  common  for  each  column  and  listed  at  the  top. 
Strain  direction  common  for  each  row  and  listed  to  on  the  left. 
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Figure  12.  Relative  potential  energy  per  mol  of  molecules  for  a-RDX  for  the 
three  cases  of  uniaxial  strain. 

Each  component  of  the  SB  potential  energy  function  is  explored  individually  to  determine  which 
energy  components  are  activated  by  each  of  the  uniaxial  strains  in  figures  13  and  14.  The  phase 
transition  from  these  plots  is  apparent  by  the  sharp  discontinuity  in  the  curves.  The  phase 
transition  induced  by  E22  leads  to  a  large  increase  in  electrostatic  energy,  shown  in  figure  13(a). 
The  large  change  in  nonbonded  energy  for  the  phase  transition  occurring  for  E22  approximately 
balance  each  other  out,  and  the  potential  energy  in  figure  12  does  not  show  a  large  break  in 
energy.  It  is  interesting  to  note  that  in  figure  14,  the  change  in  angle  and  dihedral  energy  for  E22 
is  initially  opposite  to  that  of  both  En  and  E33.  The  bond  energy  was  not  plotted  because  it  is  a 
stiff  harmonic  potential  energy  function  and  does  not  change  noticeably  across  the  phase 
transition. 


Figure  13.  (a)  Electrostatic  and  (b)  van  der  Waals  energy  per  mol  of  RDX  molecules  for  a-RDX. 
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Figure  14.  (a)  Angle  and  (b)  dihedral  energy  per  mol  of  RDX  molecules  for  a-RDX. 

3.3  Stacking  Mismatch  Simulations 

One  of  the  main  goals  of  this  research  is  to  develop  a  continuum  model  of  the  deformation 
features  in  RDX  and  other  energetic  crystals  to  determine  the  conditions  favorable  to  dislocation 
emission  and  brittle  fracture.  Dislocation  emission  leads  to  plastic  behavior,  and  brittle  fracture 
leads  to  the  creation  of  new  surface  area.  Both  of  these  mechanisms  affect  the  subsequent 
response  of  the  material  to  continued  loading.  Dislocation  emission  provides  a  means  for 
dislocations  to  build  up  at  obstacles,  as  proposed  by  Armstrong  (3).  The  obstacles  may  be  the 
result  of  limited  available  slip  systems  proposed  by  Ramos  et  al.  (5).  The  dislocations  could  also 
get  tangled  up  in  the  complex  structure  of  the  molecules,  as  in  the  steric  hindrance  model  of  Dick 
and  Ritchie  ( 4 ).  All  of  these  theories  agree  on  the  ability  of  the  material  to  emit  dislocations,  so 
it  is  important  to  understand  theoretically  under  what  conditions  RDX  will  emit  a  dislocation  and 
under  what  conditions  it  will  fracture. 

Brittle  versus  ductile  fracture  is  studied  using  energy  release  rates  of  cleavage,  Gcieav,  and 
dislocation  emission,  Gdisi-  The  energy  release  rate  for  brittle  fracture  is  understood  in  terms  of 
the  Griffith  criterion,  if  Gcieav  >  2ys,  fracture  occurs  where  ys  is  the  energy  of  a  free  surface,  and 
the  factor  of  2  comes  from  the  creation  of  two  surfaces  (20).  Ductile  behavior  is  more  complex, 
requiring  a  dislocation  to  nucleate  at  a  crack  tip  and  then  be  driven  away  from  the  crack  tip.  A 
material  behaves  in  a  ductile  manner  if  Gdis^Gdeav- 

Rice  (2)  developed  a  method  of  determining  Gdisi  by  using  the  Pierls  concept  of  creating  a  shear 
stress-atomic  displacement  interplanar  potential  that  could  be  applied  as  a  boundary  condition  to 
the  elasticity  solution  for  a  traction  free  crack  loaded  in  mode  II.  Rice  determined  Gdisi  to  be 
equal  to  a  new  solid  state  parameter  called  the  unstable  stacking  fault  energy  yus.  Within  this 
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model,  the  stress  intensity  factor,  Kn,  is  found  to  be  proportional  to  ^Yus-  This  parameter  cannot 
be  directly  determined  experimentally  but  is  readily  available  at  the  saddle  points  of  the  energy 
contours  shown  in  figures  17  and  19.  Rice  also  developed  the  theory  to  deal  with  partial  and 
mixed  mode  dislocations,  both  of  which  will  be  important  for  the  RDX  crystal.  Rice’s  model  is 
shown  to  provide  accurate  estimates  for  fee  crystals  and  metals  loaded  in  pure  mode  II.  An 
approximate  expression  is  provided  for  mixed  mode  loading. 

Analysis  using  Rice’s  model  will  incorporate  atomistic  effects  into  the  continuum  description  of 
the  material’s  deformation.  This  will  provide  the  ability  to  predict  the  deformation  mechanism 
in  RDX  based  on  a  variety  of  loading  conditions.  Rice’s  model  can  also  incorporate  temperature 
dependence  to  predict  the  temperature  at  which  the  deformation  mechanism  under  an  applied 
mode  of  loading  switches  from  brittle  to  ductile. 

The  SB  potential  will  be  used  to  determine  the  generalized  stacking-mismatch  energy  contours, 
or  y-surface,  for  the  crystallographic  planes  of  a-RDX.  The  y-surface  refers  to  an  energy  contour 
and  is  not  related  to  the  y-RDX  phase.  The  planes  to  be  studied  are  the  experimentally 
determined  slip  planes  (5).  The  y-surface  will  be  used  to  determine  parameters  for  continuum 
level  models  of  fracture  and  dislocation  emission,  such  as  Rice’s  dislocation  nucleation  model 
(2).  This  section  will  provide  details  of  the  stacking  fault  simulations  that  have  been  done  for  the 
(010)  and  (001)  planes. 

Figure  15(a)  shows  the  initial  crystal  lattice  used  to  create  stacking  faults  on  the  (010)  slip  plane. 
The  lattice  shown  contains  3x30x3  unit  cells  of  a-RDX.  This  lattice  will  be  used  to  create  a 
stacking  fault  on  the  (010)  plane  by  displacing  the  top  half  of  the  lattice  by  a  vector  u=ua+uc 
relative  to  the  bottom  half  along  the  (010)  slice  plane  shown  in  figure  15(b).  The  purpose  of  the 
simulations  is  to  create  a  single  infinite  stacking  mismatch  across  the  slip  plane  surrounded  by 
the  bulk  material.  The  infinite  slip  plane  is  created  using  periodic  orthorhombic  boundary 
conditions  in  the  a-  and  c-directions.  Creating  a  single  stacking  mismatch  across  the  slip  plane 
surrounded  by  bulk  material  in  the  b-direction  is  more  difficult  and  requires  the  creation  of  a 
sandwiched  structure  of  active  molecules  surrounded  by  frozen  molecules  to  simulate  the  bulk. 

A  vacuum  is  created  at  the  top  of  the  structure  so  as  to  not  create  a  second  stacking  fault  due  to 
periodic  boundary  condition  in  the  b-direction.  The  frozen  atoms  shield  the  effect  of  the  free 
surface  from  the  active  molecules  at  the  slip  plane  and  also  constrain  the  motion  of  active 
molecules.  A  similar  sandwich  structure  and  simulation  procedure  was  used  by  Bedrov  et  al.  to 
study  the  development  of  y-RDX  in  the  vicinity  of  shock  loads  (10). 

The  stacking  mismatches  for  the  (010)  plane  are  created  by  the  following  procedure: 

•  Anneal  3x3x3  lattice  using  NPT  ensemble  at  P=0  and  T=25K.  NPT  ensemble  with 

orthorhombic  boundary  conditions  is  used  to  constrain  lattice  vectors  to  remain  orthogonal 
during  simulation.  Annealing  process  involves  a  series  of  25K  thermalizations  followed  by 
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conjugate  gradient  force  minimizations.  The  minimum  energy  configuration  from  the 
annealing  process  is  used  in  the  remainder  of  steps. 

•  Anneal  minimized  configuration  from  step  1  using  an  NVT  ensemble  and  keep  minimum 
energy  configuration. 

•  Use  the  averaged  3x3x3  system  to  construct  a  3x30x3  crystal  lattice  by  repeating  it  10 
times  in  the  b-direction. 

•  In  the  new  3x30x3  superlattice,  freeze  the  atomic  positions  of  the  top  14  and  bottom  14  unit 
cells  in  the  b-direction 

•  Add  a  4  nm  vacuum  region  to  the  top  of  the  simulation  cell  by  increasing  the  b-lattice  by 
4  nm. 

•  Apply  strain  in  the  b-direction  to  the  layer  of  flexible  molecules  on  top  of  and  below  the 
stacking  fault  to  displace  them  away  from  each  other  along  the  stacking  fault.  This  is  done 
to  keep  the  atomic  coordinates  of  the  molecules  from  coinciding  at  the  start  of  the 
simulation  from  the  displacements  applied  in  the  next  step. 

•  Displace  the  top  half  of  molecules  by  a  fraction  of  the  |a|  and  |c|  lattice  lengths  so  as  to 
create  a  grid  of  100  displaced  configurations  on  the  (010)  plane. 

•  Run  NVT  annealing  simulations  at  T=25  K  for  each  new  stacking  mismatch  configuration. 

•  For  each  simulation,  the  minimum  potential  energy  provides  a  single  scalar  energy  value 
on  the  10x10  y-surface  plotted  in  figure  17. 

•  The  minimized  atomic  coordinates  of  the  flexible  layers  are  used  to  find  the  structure  and 
molecule  conformation  of  the  stacking  mismatches. 

Note:  The  same  procedure  is  used  to  create  the  (001)  stacking  mismatches,  except  that  unit  cells 
are  stacked  in  the  (001)  direction. 
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Figure  15.  (a)  3x30x3  a-RDX  sandwich  structure  where  the  blue  region  is  4-nm  vacuum,  red  are  frozen  atoms,  and 
green  are  frilly  flexible  molecules  with  slice  plane  through  center  of  green,  (b)  Schematic  of 
displacement,  u  =  ua  +  uc,  of  upper  lattice  on  glide  plane. 


The  cut  plane  for  the  (010)  stacking  fault  is  shown  in  figure  16.  For  displacement  in  the  a- 
direction,  no  molecules  come  into  contact  with  each  other.  Displacement  in  the  c-direction 
requires  the  molecules  across  the  slip  plane  to  move  in  the  b-direction  in  order  to  move  out  of 
each  other’s  way.  This  strains  the  lattice  locally,  and  potential  energies  for  these  stacking  faults 
are  large.  This  steric  effect  of  large  molecules  blocking  available  slip  systems  is  common  to 
energetic  crystals  and  molecular  crystals  in  general. 


Figure  16.  Side  view  of  the  stacking  fault  cut  plane  in  blue  for  the  (010)  slip  plane. 
Molecules  above  the  blue  line  are  moved  in  the  c-  and  a-directions. 
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The  y-surface  is  shown  in  figure  17(a)  for  the  (010)  slip  plane  shown  in  figure  16.  Large  energy 
barriers  are  shown  by  the  red  contours  for  displacement  in  the  c-direction,  also  shown  by  the  red 
line  in  the  energy  plot  in  figure  17(b).  A  much  smaller  energy  barrier  is  encountered  for 
displacements  in  the  a-direction,  shown  by  the  green  line  in  figure  17(b).  The  dashed  line  in 
figure  17(b)  is  the  free  surface  energy  for  the  (010)  plane  multiplied  by  two. 


Figure  17.  (Left)  y-surface  for  (010)  stacking  fault,  and  (right)  slices  of  y-surface  along  a-  (green)  and  c-axis 
(red).  Saddle  point  and  free  surface  energies  labeled 


The  cut  plane  for  the  (001)  stacking  fault  is  shown  in  figure  18.  Again,  for  displacement  in  the 
a-direction,  no  molecules  come  into  contact  with  each  other.  Displacement  in  the  b-direction 
requires  the  molecules  across  the  slip  plane  to  move  in  the  c-direction  in  order  to  move  out  of 
each  other’s  way.  This  strains  the  lattice  locally,  and  potential  energies  for  these  stacking  faults 
are  large. 


Figure  18.  Side  view  of  the  stacking  fault  cut  plane  in  green  for  the  (001) 

stacking  fault.  Molecules  above  the  right  of  the  line  are  moved  in 
the  b-  and  a-directions. 
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The  y-surface  is  shown  in  figure  19(a)  for  the  (010)  slip  plane  shown  in  figure  18.  Smaller 
energy  barriers  occur  for  the  sterically  hindered  displacement  in  the  b-direction  on  this  slip  plane 
than  for  the  (001)[00 1]  slip  system  shown  by  the  red  line  in  figure  17(a).  A  much  smaller  energy 
barrier  is  encountered  for  displacements  in  the  a-direction,  shown  by  the  green  line  in  figure 
19(b).  The  dashed  line  in  figure  19(b)  is  the  free  surface  energy  for  the  (001)  plane  multiplied  by 
two.  Results  for  the  (001)  and  (010)  slip  planes  are  summarized  in  table  5. 


Figure  19.  (Left)  y-surface  for  (001)  stacking  fault,  and  (right)  slices  of  y-surface  along  a-  (green)  and  b-axis 
(red).  Saddle  point  and  free  surface  energies  labeled 


Table  5.  Summary  of  slip  systems  and  prediction  of  deformation  mechanism  based  on  Rice’s  analysis  (2) 


Summary  of  Slip  Systems 

Slip  Systems 

yus  (mJ/m2) 

2ys  (mJ/m2) 

ratio 

This  Work 

Experimental'1 

(010)[001] 

564 

404 

1.40 

Brittle 

Cross-slip 

(010)[100] 

215 

404 

0.53 

Ductile 

Slip  System 

(001)[010] 

398 

427 

0.93 

Ductile 

Cleavage  Plane 

(001)[100] 

250 

427 

0.59 

Ductile 

Cleavage  Plane 

aRamos  et  al.  (5) 

Table  5  also  presents  predictions  based  on  Rice’s  analysis  to  determine  the  probable  mode  of 
deformation  in  pure  mode  II  loading  by  comparing  yus  to  2ys.  Based  on  the  y-surfaces  found  in 
this  work,  only  the  (0 1 0)[00 1  ]  slip  system  is  found  to  be  brittle.  Ramos  et  al.  (5)  find  this  to  be  a 
cross-slip  plane,  meaning  that  it  is  experimentally  ductile.  The  ductility  Ramos  et  al.  observe  for 
this  plane  from  nano-indentation  could  be  due  to  cross  slip  after  initial  displacement  to  the  stable 
stacking  fault  configuration  at  (0 1 0)[002]  occurs.  The  (0 1 0)[00 1  ]  system  may  also  not  evolve 
directly  in  the  [001]  direction  and  may  vary  slightly  as  to  avoid  the  large  energy  barrier  at  [002], 
shown  by  the  red  contour  in  figure  17(a). 
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Ramos  et  al.  also  reported  the  (001)  plane  to  be  a  cleavage  plane,  but  in  this  work  it  is  shown  to 
be  ductile.  Ramos  et  al.  cite  other  micro-indentation  hardness  studies  as  their  source  for 
cleavage  on  this  plane.  Microindentations  studies  on  RDX  crystals  probably  apply  mixed  mode 
loading,  which  would  aid  in  the  development  of  fracture.  The  ratio  of  yus  to  2ys  on  this  plane  is 
also  slightly  larger  than  that  found  for  the  (010)  plane,  and  it  should  be  slightly  more  susceptible 
to  brittle  fracture. 

Other  experimental  slip  systems  such  as  (02 1)[100]  and  (01 1)[  1 00]  reported  by  Ramos  et  al. 
need  to  be  constructed  and  studied.  One  difficulty  in  these  simulations  is  correctly  determining 
the  slip  plane.  The  most  probable  slip  plane  for  the  (021)[100]  slip  system  is  shown  in  figure 
20(a)  by  the  box  contained  between  the  blue  and  green  dashed  lines.  The  dashed  lines  in  figure 
20  represent  the  slip  planes.  However,  determining  where  these  slip  planes  intersect  unit  cells, 
and  which  molecules  lie  on  which  side  of  the  slip  system,  is  difficult.  This  is  the  case  for  the 
(01 1)[  1 00]  slip  system,  and  one  possible  unit  cell  is  shown  between  the  dashed  lines  in  figure 
20(b).  All  of  the  unit  slip  vectors  activated  experimentally  by  Ramos  et  al.  are  in  the  [100] 
direction,  which  is  not  surprising  because  this  is  the  only  Burgers  vector  for  which  slip  can  occur 
without  the  molecules  running  into  one  another — i.e.,  this  is  the  least  sterically  hindered  Burgers 
vector.  Geometrically,  dislocations  favor  the  smallest  lattice  dimension  because  the  dislocation 
energy  scales  by  |b|  ,  where  |b|  is  the  length  of  the  Burgers  vector. 


Figure  20.  Other  experimental  a-RDX  slip  systems  and  with  dashed  representing  slip  planes. 


4.  Summary  and  Conclusions 


The  high  pressures  that  occur  in  detonation  fronts  interacting  with  energetic  crystals  are  in  the 
pressure  range  where  phase  transitions  may  occur.  The  properties  of  these  new  solid  phases  are 
unlike  those  found  at  ambient  conditions  and  may  play  a  vital  role  in  the  development  of  hot 
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spots.  This  work  used  pressurization  simulations  to  show  the  ability  of  the  SB  potential  to 
predict  both  the  high  and  low  pressure  (a  and  y)  phases  of  RDX  and  transition  from  the  high 
pressure  phase  to  the  low  pressures  phase  (y-  to  a-RDX).  The  reverse  phase  transition  was  not 
observed  using  NPT  simulations,  and  so  NVT  simulations  of  uniaxial  strain  were  conducted  for 
the  low  pressure  a-RDX  phase  to  determine  if  the  phase  transition  was  dependent  on  the  load 
orientation.  Applying  uniaxial  strain  along  the  c-axis  caused  the  material  to  convert  from  the  a- 
to  y-RDX.  This  crystal  dependence  to  load  orientation  is  of  interest  because  experimental  data 
by  Cawkwell  et  al.  (9)  show  the  orientation  dependence  of  shock  properties,  and  their  simulation 
work  determined  that  the  dislocation  structures  involved  the  phase  transition  from  a-  to  y-RDX. 

The  deformation  mechanisms  of  RDX  are  of  interest  to  the  energetic  materials  community 
because  of  their  role  as  a  localization  feature  in  perfect  crystals  under  low  strain  rate  loading, 
which  occurs  during  accidental  mishandling.  This  work  provides  an  initial  estimate  and  method 
of  predicting  if  dislocation  nucleation  or  brittle  fracture  will  occur  based  on  the  continuum  level 
stress  intensity  factors  using  Rice’s  model.  These  atomistic  level  details  of  the  deformation 
mechanism  that  are  active  for  low  strain  rate  loads  are  lacking  in  both  simulation  and 
experimental  literature.  Cawkwell  et  al.  (9)  report  that  according  to  the  steric  hindrance  model, 
deformation  along  orientations  containing  active  or  inactive  slip  systems  are  either  insensitive  or 
sensitive,  respectively.  By  determining  the  y-surface  for  the  different  experimental  and 
crystallographic  planes,  as  was  done  in  this  work,  Rice’s  analysis  can  be  used  to  determine  the 
active  and  inactive  slip  systems,  and  from  that,  determine  orientations  to  which  the  material  is 
most  sensitive  to  shock. 

The  future  work  will  provide  a  methodology  for  studying  crack  tip  dislocation  emission  and 
brittle  fracture  of  RDX,  which  is  also  applicable  to  other  molecular  crystalline  materials  held 
together  by  nonbonded  forces.  The  future  work  will  be  useful  to  the  modeling  community  by 
providing  information  on  the  deformation  features  that  need  to  be  captured  in  the  continuum 
level  hydrocodes.  These  hydrocodes  are  important  to  the  development  and  testing  of  new 
military  munitions  to  determine  their  effectiveness  and  safety.  By  better  understanding  the 
deformation  and  fracture  process,  the  manufacturing  process  of  munitions  can  be  made  safer  and 
provide  more  effective  munitions.  Explosive  fills  are  sometimes  compressed  or  machined  for 
munitions,  and  defects  introduced  during  this  stage  may  affect  both  their  performance  and  safety. 
This  is  also  true  in  the  pharmaceutical  industry,  where  several  drugs  of  interest  use  molecular 
crystals  as  active  ingredients.  The  molecular  crystals  are  pressed  into  pill  form,  and  fractures 
introduced  during  this  step  will  affect  the  rate  at  which  the  body  is  able  to  break  them  down  and 
absorb  them. 
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Abstract 


A  domain-specific  term  base  may  be  useful  not  only  as  a  resource  for  written  and  oral 
translation,  but  also  for  Natural  Language  Processing  (NLP)  applications,  text  retrieval, 
document  indexing,  and  other  knowledge  management  tasks.  The  objective  of  this  investigation 
was  to  explore  the  use  of  alternative  terminology  extraction  methods  to  refine  and  validate  an 
existing  military-specific  bilingual  dictionary.  A  series  of  semi-automatic  methods  was 
implemented  to  distill  the  existing  term  list  by  removing  redundancies,  resolving  spelling 
variations,  and  separating  individual  expressions.  Once  the  internal  clean-up  was  completed,  we 
compared  two  methods  drawn  from  the  terminology  extraction  literature  in  order  to  validate 
terms  as  military-specific  and  to  propose  a  candidate  list  of  non-specific  terms  for  exclusion — 
term  frequency  calculations  and  terminology  extraction  lists.  In  this  investigation,  we  wanted  to 
find  the  best  procedure  to  extract  domain-specific  terms  for  a  low-resource  domain;  to 
demonstrate  that  terminology  extraction  methods  can  be  used  to  validate  and  refine  a  domain- 
specific  dictionary;  and  to  provide  the  final,  refined  dictionary  as  a  term  base  to  support 
customization  of  machine  translation  systems  for  the  military  domain. 
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1.  Introduction 


Especially  since  the  2001  entrance  of  the  United  States  into  the  war  in  Afghanistan,  foreign 
language  translation  has  become  increasingly  necessary  yet  still  is  not  sufficiently  resourced. 
Although  human  translators  often  provide  high-quality  work,  that  work  can  be  costly  and  time 
consuming  given  that  it  is  difficult  to  find  qualified  bilingual  language  experts  across  all  needed 
domains.  This  lack  of  quick  translation  along  with  advances  in  the  information  technology  field 
has  prompted  research  into  and  use  of  semi-automatic  machine  translation  (MT)  methods  to 
support  human  translators.  Whereas  word-to-word  translation  in  specialized  domains  may  be 
straightforward  (e.g.,  stethoscop e-estetocopio)  given  a  language  expert  or  a  bilingual  dictionary, 
the  difficulty  lies  with  multi-word  expressions — with  recognizing  phrases  that  are  in  fact 
technical  terms  (“field  of  fire”)  and  need  to  be  treated  as  entities,  and  with  finding  their 
counterparts  in  the  other  language,  where  the  phrase  may  or  may  not  have  the  equivalent  number 
of  words. 

Over  the  last  10  years,  tools  to  enable  automatic  extraction  of  term  bases  have  been  developed, 
which  speed  the  process  of  deriving  term  bases  from  a  collection  of  documents  in  a  domain  of 
interest.  A  domain-specific  term  base  may  be  useful  not  only  as  a  resource  for  written  and  oral 
translation,  but  also  for  Natural  Language  Processing  (NLP)  applications,  text  retrieval  (7), 
document  indexing,  and  other  knowledge  management  tasks.  The  National  Virtual  Translation 
Center  (NVTC),  an  organization  under  the  Lederal  Bureau  of  Investigation,  was  established  in 
Lebruary  2003  for  the  exact  purpose  of  “providing  timely  and  accurate  translations  of  foreign 
intelligence  for  all  elements  of  the  intelligence  community  (2).”  In  September  of  that  year,  an 
electronic  compilation  of  8953  terms  with  their  translations  was  published  by  M.  Green  for  the 
NVTC,  under  the  title  Iraqi  Military  English-Arabic  Arabic-English  Dictionary.  While  the 
sources  of  these  translated  terms  and  the  purpose  of  the  dictionary  are  unclear,  it  has  been  used 
successfully  to  support  improved  MT. 


2.  Examining  the  NVTC  Bilingual  Military  Dictionary 


Searching  through  the  original  term  list,  we  found  many  internal  discrepancies  and 
inconsistencies  that  suggested  that  the  term  base  may  have  been  developed  by  several  authors 
and  provided  rapidly  to  the  field  for  urgent  needs  without  opportunity  for  quality  assurance. 
These  internal  issues  would  pose  problems  with  its  use  in  computational  linguistics.  The 
problems  include  the  following: 
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1 .  Alignment  and  spacing  errors 

a.  White  space  preceding  the  expression  alters  its  place  when  ordered  alphabetically. 

b.  White  space  trailing  the  expression  can  introduce  two  entries  from  the  same  expression: 

i.  Example:  One  entry  would  be  given  as  “Flank”  while  the  other  would  be  provided 
as  “Flank”  and  they  would  have  the  same  Arabic  translation. 

2.  Thirty-three  duplicate  entries 

a.  These  entries  are  exactly  the  same  in  both  Arabic  and  English;  therefore,  the  duplicates 
can  be  removed. 

3.  Three  variations  of  the  same  word 

a.  The  dictionary  would  include  two  non-identical  English  entries  with  the  identical 
Arabic  translation: 

i.  Example:  “Light  antiaircraft”  and  “Light  anti-aircraft”  had  the  same  Arabic 
translation  “AA  ^iyu=,  yio”  and  “AA  ^j!yu=,  A”. 

b.  For  the  purposes  of  this  project,  both  entries  were  used,  but  at  the  end  of  the 
investigation,  only  the  most  commonly  used,  grammatically  correct  entry  was  included 
in  the  dictionary. 

4.  Five  misspellings 

a.  Example:  “Airconditioned  shelter”  should  be  “Air-conditioned  shelter”. 

b.  When  air-conditioned  is  listed  as  its  own  entry,  it  has  the  appropriate  spelling,  but  when 
combined  with  another  word,  it  is  spelled  incorrectly. 

5.  An  unnecessary  symbol,  o,  was  included  after  three  English  entries. 

6.  For  computational  linguistic  purposes,  tokenizations  would  have  to  be  performed  on  the 
following  collections:  parentheses  (622),  ampersands  (15),  and  slashes  (166).  A  blank 
space  was  inserted  where  the  original  character  was  located. 

Arabic  experts  looked  at  a  random  sample  of  the  existing  terminology  that  I  proposed  as 
representative  and  noted  that  (1)  the  terminologies  were  of  many  cultural  dialects,  but  mainly 
Standard  Modem  Arabic,  and  (2)  the  Arabic  translation  of  general  English  words  did  not  have  a 
military-specific  connotation,  suggesting  that  the  term  does  not  belong  in  the  dictionary.  Since 
we  are  simply  focusing  on  the  English  portion  of  the  term  base,  its  bilingual  nature  does  not 
really  enter  into  the  processes  used  to  refine  the  dictionary  at  this  time.  Further  research  is 
needed  for  the  Iraqi-Arabic  portion. 
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3.  Internal  Clean-Up 


In  order  to  make  the  existing  term  base  ready  for  computer  intervention,  several  changes  had  to 
be  made  (noted  in  figure  1).  Using  a  Perl  script,  we  found  that  the  original  NVTC  term  base  had 
8953  entries  with  the  following  breakdown: 


WPL:  1 

AOL:  1832 

WPL:  9 

AOL: 

WPL:  2 

AOL:  4795 

WPL:  10 

AOL: 

WPL:  3 

AOL:  1591 

WPL:  1 1 

AOL: 

WPL:  4 

AOL:  440 

WPL:  13 

AOL: 

WPL:  5 

AOL:  182 

WPL:  14 

AOL: 

WPL:  6 

AOL:  83 

WPL:  16 

AOL: 

WPL:  7 

AOL:  24 

WPL:  18 

AOL: 

WPL:  8 

AOL:  15 

WPL:  19 

AOL: 

WPL:  Words  per  Line 
AOL:  Amount  of  Lines 


Figure  1.  Internal  correction  process. 


Once  we  became  familiar  with  the  term  base,  we  determined  that  it  had  to  be  altered  in  order  to 
accurately  process  the  material.  The  list  of  problems  identified  in  the  introduction  was  used  to 
refine  existing  text.  First,  the  terms  were  alphabetized.  Entries  that  had  unnecessary  preceding 
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white  space  were  fixed.  Microsoft  Office  was  unable  to  remove  trailing  white  spaces,  so  Perl  was 
used  for  this  purpose.  The  code  removed  all  white  space  after  each  string  in  the  text  file  and 
replaced  the  new  entry  in  the  dictionary. 

Once  the  alignment  and  spacing  errors  were  corrected,  both  a  Perl  script  and  Conditional 
Formatting  within  Microsoft  Excel  were  used  to  identify  all  exact  matches  within  the  column  of 
terms.  Both  methods  identify  a  total  of  33 1  duplicates  in  the  English  portion.  Taking  the  entire 
dictionary  into  context,  there  were  33  duplicate  entries  (some  entries  were  found  three  separate 
times);  therefore,  37  entries  were  removed. 

In  response  to  the  variations  among  words  in  the  dictionary,  we  decided  to  include  both  entries  to 
find  the  most  common  spelling  in  order  to  eliminate  one  of  the  entries  later  in  the  project. 
Misspellings  were  then  corrected  to  help  reinforce  standardization  of  the  term  base.  We  also 
removed  the  unnecessary  symbol  following  three  of  the  entries. 

Entries  with  two  separate  terms  combined  and  submitted  as  one  entry  were  noted  (i.e., 
antiaircraft/artillery,  director/directorate).  These  submissions  should  be  separated  into  two  entries 
for  the  purpose  of  accessibility  in  the  field,  and  in  our  term  frequency  method,  exact  string 
matching  is  essential  for  accurate  results.  Therefore,  all  entries  with  gratuitous  explanations  and 
definitions  following  the  term  were  removed.  A  Microsoft  Excel  macro  was  employed  to 
eliminate  all  items  within  parentheses. 

Once  these  alterations  were  completed,  the  new  term  base  consisted  of  the  following  breakdown: 


WPL:  1 

AOL:  1832 

WPL:  9 

AOL: 

WPL:  2 

AOL:  4795 

WPL:  10 

AOL: 

WPL:  3 

AOL:  1591 

WPL:  1 1 

AOL: 

WPL:  4 

AOL:  440 

WPL:  13 

AOL: 

WPL:  5 

AOL:  182 

WPL:  14 

AOL: 

WPL:  6 

AOL:  83 

WPL:  16 

AOL: 

WPL:  7 

AOL:  24 

WPL:  18 

AOL: 

WPL:  8 

AOL:  15 

WPL:  19 

AOL: 

WPL:  Words  per  Line 
AOL:  Amount  of  Lines 
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4.  Method  One:  Frequency  Count 


This  proposed  method  to  collection  a  set  of  domain-specific  terminology  is  based  on  the 
principle  of  Term  Frequency-Inverse  Document  Frequency  (TF-IDF).  As  tested  in  An 
Unsupervised  Approach  to  Domain-Specific  Term  Extraction  (3),  the  principle  behind  frequency 
counting  is  the  idea  that  certain  terminology  will  generally  occur  with  a  higher  frequency  within 
domain-specific  documents  as  opposed  to  in  a  general  corpus.  This  theory,  however,  has  its 
limitations.  Single  word  terminology  is  much  more  difficult  to  access  based  on  the  occurrences 
of  homographs.  In  the  NVTC’s  dictionary  for  example,  the  entry  “brief’  could  be  found  in 
several  different  contexts.  In  a  military  sense,  the  term  can  be  used  as  a  verb  to  summarize  or 
give  preparatory  information  to  Soldiers,  but  in  a  general  connotation,  it  could  be  used  as  an 
adjective  or  noun  to  describe  duration  and  length. 

4.1  Input 

A  domain-specific  corpus  of  2619  documents  was  then  created  by  collecting  various  military 
documents  from  a  variety  of  sources.  The  documents  selected  were  chosen  because  of  their 
translated  nature;  if  a  document  was  important  enough  to  military  use  that  it  was  translated  into 
Arabic,  then  its  extracted  terminology  is  most  likely  vital  to  a  bilingual  dictionary.  Thirteen  items 
from  the  Ranger  Handbook,  one  item  from  field  manual  3-2 1.10,  and  five  items  from  field 
manual  7-8  were  selected,  along  with  93  documents  from  the  Combating  Terrorism  Center’s 
Harmony  Database  of  Released  Documents  (CTC)  and  2507  items  from  an  Iraqi  database  from 
ARL’s  holdings.  The  CTC  at  West  Point,  dedicated  to  scholarly  research  and  policy  analysis  to 
examine  combat  terrorism,  published  a  series  of  letters,  reports,  and  al-Qa’ida-related  documents 
captured  during  the  War  on  Terror  for  public  access.  This  is  important  to  our  corpus  as  a  first¬ 
hand  account  of  events  in  Afghanistan,  elucidating  al-Qa’ida’s  actions  and  weaknesses.  The  Iraqi 
training  material  consists  of  PowerPoint  training  materials,  scripts,  and  guides  to  a  variety  of 
field  situations. 

4.2  Output 

The  goal  of  this  method  was  to  take  the  internally  cleaned  dictionary  and  use  exact  string 
matching  to  search  through  the  corpus  for  the  number  of  occurrences  of  each  term.  Because  of 
the  extensive  nature  of  the  corpus,  we  used  a  Hadoop  cluster,  a  programming  framework 
designed  for  large-scale  computational  use,  to  expedite  the  process.  Before  processing  the  data, 
all  the  documents  (Acrobat  Reader,  Microsoft  Word,  Microsoft  Excel,  and  Microsoft 
PowerPoint)  were  converted  into  text  files  with  the  help  of  an  online  converter.  The  Iraqi  training 
documents  could  not  be  easily  converted,  however,  because  of  the  high  number  of  subfolders 
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within  each  main  folder.  Again  using  Perl,  we  renamed  all  documents,  changing  spaces  to  dashes 
and  ampersands  to  underscores,  and  moved  all  documents  to  one  large  folder,  which  helped  ease 
the  conversion  of  the  files. 

Once  all  target  files  were  converted,  they  were  processed  with  the  servers  searching  for  exact 
string  matches  based  on  the  dictionary’s  terms.  The  process  resulted  in  two  Excel  files 
summarizing  the  findings.  The  first,  “Word  Count”  (table  1),  was  a  list  of  all  keywords,  the 
number  of  occurrences  in  the  corpus,  and  on  average  how  many  times  that  keyword  appeared  per 
document.  The  second  file,  “Doc  Count”  (table  2),  consisted  of  a  list  of  each  document,  the 
number  of  key  words  in  the  document,  and  the  average  number  of  times  a  keyword  appeared. 


Table  1.  Word  Count  chart  excerpt. 


Term 

No.  of  Times  Term 
Appears  in  Corpus 

Map  reconnaissance 

16 

Fallout 

16 

Psychological  warfare 

16 

Stud 

16 

Barrel  assembly 

16 

Medical  unit 

15 

Table  2.  Doc  Count  chart  excerpt. 


Document 

No.  of  Terms  in  Dictionary  that 
Appear  in  Corpus 

Iraqi-Training-Disk_S3_MOUT_  Infantry-Rifleman- 
Course-Handout-Booklet-2003.txt 

462 

Iraqi-Training-Disk  ca-documents  instant-lessons-of- 
iraq-war.txt 

458 

AFGP-2002-600092-Trans-Meta.txt 

448 

Iraqi-Training-Disk  ca-documents  SASO-handbook.txt 

434 

AFGP-2002-600088-Trans-Meta.txt 

371 

AFGP-2002-600053-Trans-Meta.txt 

361 
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The  results  from  the  TF  method  indicated  that  the  most  common  terms  were  as  follows: 


One-word  entries 

Enemy 

8622 

Support 

5254 

Commander 

4889 

Operations 

4874 

Two-word  entries 

First  aid 

448 

Armed  forces 

389 

Indirect  lire 

340 

Warning  order 

316 

Three-word  entries 

Course  of  action 

364 

Command  and  control 

310 

Chain  of  command 

306 

Concept  of  operations 

216 

The  results  support  Zipf  s  Law  ( 4 )  that  term  length  is  inversely  proportional  to  its  number  of 
occurrences  in  a  corpus.  Zipf  s  Law  will  become  an  important  factor  in  the  term  extraction 
process.  We  found  29.68%  of  all  terms  in  the  dictionary  with  a  frequency  of  one  or  more  in  the 
corpus  and  26. 13%  of  those  appeared  more  than  once. 


5.  Method  Two:  Terminology  Extraction 


The  goal  of  terminology  mining  or  extraction  is  to  collect  a  list  of  domain-pertinent  terms  from  a 
given  corpus.  For  the  purposes  of  this  investigation,  the  online  extraction  tool  TermExtractor  (5), 
developed  by  the  Linguistic  Computing  Laboratory  of  the  University  of  Roma,  was  used  to 
determine  what  percentage  of  the  extracted  term  list  overlapped  with  the  existing  military  bank. 

The  terms  that  appear  in  both  corpora  are  then  added  to  a  proposed  list  of  confirmed  dictionary 
entries.  Figure  2  shows  the  TermExtractor  pipeline. 


51 


Figure  2.  TermExtractor  pipeline  (5). 

To  ensure  consistency  in  our  results,  we  used  the  same  corpus  as  a  reference  throughout  the 
entire  project.  We  submitted  the  same  corpus  of  2619  documents  as  in  the  TF  method  to  be 
processed  for  specificity.  TermExtractor  uses  input  documentation  to  extract  statistically  relevant 
terminology  through  the  use  of  chuncking  and  document  parsing,  as  well  as  by  filtering 
unecessary  information.  These  filters  eliminate  stopwords  such  as  “the,  as,  is,  for”  and  general 
terminology  that  does  not  indicate  domain-specificity.  The  extraction  tool  filters  non- 
terminological  strings  through  its  evalution  of  the  following: 

•  Domain  Pertinence:  High  (numerical  value)  means  a  term  is  frequent  in  the  domain  of 
interest  and  is  much  less  frequent  in  the  other  domains  used  for  contrast  ( 6 ): 

DRDi  (t)  =  -  X  P~  (t  /dk)  log  (P~  (t  /dk ))  =  X  norm_freq(t,dk)  log  (norm  _freq(t,dk)) 

•  Lexical  Cohesion:  The  degree  to  which  the  terms  adhere  to  one  another  within  a  string. 
This  proved  more  effective  than  other  measures  of  cohesion  (6).  The  resulting  numerical 
value  is  high  if  the  words  within  a  string  occur  more  often  with  one  another  rather  than 
alone  in  a  corpus.  The  minimum  was  set  to  0.05. 

•  Structural  Relevance:  When  a  title  or  subtitle  is  composed  of  domain-specific  terms,  then 
its  importance  is  increased  by  some  factor  x.  Highlighted,  bolded,  and  italicized  items  are 
also  included  (x=5  for  highlighted,  capitalized,  underlined,  colored,  smallcaps,  italicized, 
and  bolded  terms,  and x=10  for  titles  and  abstract  content). 

•  Miscellaneous:  A  set  of  heuristics  are  applied  to  increase  computational  performance  by 
removing  generic  articles  and  terminology,  detecting  misspellings,  distinguishing  part  of 
speech,  extracting  uni  gram  terminology,  and  detecting  abbreviations. 

The  extraction  tool  also  sets  up  contrastive  corpora  to  eliminate  common  terminology  that  may 
be  relevant  to  the  specific  domain  but  not  entirely  of  that  domain.  These  corpora  include  the 
following: 

•  Brown  Corpus  (3634  terms) 

•  Medicine  (228 1  terms) 

•  Computer  Networks  (16335  terms) 

•  Sports  (1020  terms) 
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•  Tourism  (55590  terms) 

•  Wall  Street  Journal — Economy  (3606  terms) 

Although  these  terminology  banks  are  not  specifically  indentified,  it  is  important  to  set  up  some 
contrasting  corpora  to  eliminate  general  terminology  and  possibly  create  a  proposed  list  of  terms 
for  expulsion. 

5.1  First  Investigation 

In  the  first  investigation,  the  corpus  was  submitted  without  any  restrictive  measures  to  find  the 
percentage  of  extracted  terminology  that  would  overlap  with  the  existing  term  bank.  Given 
Zipf  s  Law  ( 4 ),  the  frequency  distribution  of  word  length  is  exponential;  this  means  that,  in 
accordance  with  a  general  corpus,  a  unigram  (one  word  term)  is  far  more  likely  to  occur  than  a 
bigram  and  a  trigram,  and  so  forth.  Due  to  time  constraints,  this  law  was  employed,  so  any  term 
that  exceeded  three  words  was  considered  domain-specific  because  of  its  exclusivity  to  a 
particular  domain.  For  all  one-  to  three-word  terms,  3605  words  occurred  in  both  the  term 
extraction  list  and  the  NVTC  dictionary.  This  indicates  that  40.27%  of  the  dictionary  is 
supported  by  this  method;  43.87%  of  all  unigrams,  bigrams,  and  trigrams. 

5.2  Second  Investigation 

For  the  second  investigation,  we  entered  the  corpus  and  entered  the  existing  term  bank  as  a 
restrictive  option.  The  extracted  terminology  from  this  trial  excludes  all  terms  in  the  dictionary  in 
its  proposed  terminology  list.  At  this  point  in  the  process,  a  human  validator  is  required  to 
identify  the  reliability  of  the  extracted  list.  I  randomly  sampled  10%  of  the  terms  (648  items)  and 
a  subject  matter  expert  evaluated  this  list,  indicating  whether  the  term  was  military-unique 
(18.06%  of  the  sample)  and  highlighting  the  spelling  errors  (24.07%).  Table  3  is  an  excerpt  of 
the  described  process,  with  its  proposed  spelling  corrections  in  column  four. 


Table  3.  Methods  comparison  to  dictionary. 


Term 

Military  Specific 

Spelling  Error 

Possible  Correction 

improvised  sling 

Yes 

include-ytank  crewmembers 

Yes 

Yes 

"including  tank  crewmembers" 

includingthe  regulationsandlaws 

Yes 

"including  the  regulations  and  laws" 

indecision  recklessness 

index  contour  line 

Yes 

This  list  will  be  used  later  as  a  basis  for  what  could  be  added  to  the  dictionary.  In  order  to  refine 
the  extracted  list  of  terms,  the  same  course  of  action  can  be  taken  as  for  the  NVTC  dictionary. 
The  possible  list  of  terms  can  be  evaluated  for  its  frequency  in  a  new  corpus  and  a  new  list  of 
terms  can  be  extracted  and  compared  for  its  similarities. 
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6.  Results 


Although  time  constraints  did  not  allow  the  full  investigation  to  be  executed,  the  original  term 
base  can  be  successfully  modified  and  refined  after  comparing  the  dictionary  with  a  general 
corpus  and  using  IDF.  The  first  portion  of  figure  3  indictates  the  overlap  between  the  orignial 
NVTC  dictionary  and  the  results  of  the  two  methods.  It  appears  that  the  TF  method  produces  a 
better  comparison  to  refining  an  existing  military  term  base,  but  the  term  extraction  method 
contributed  as  well.  The  second  portion  of  figure  3  indicates  the  overlap  between  the  TF  method 
and  the  term  extraction  method. 


Figure  3.  Comparison  to  dictionary. 

In  this  study,  27.06%  of  terms  that  appeared  with  high  frequency  also  appeared  in  the  term 
extraction  list. 
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In  addition  to  assessing  the  term  frequency  of  the  dictionary  when  paired  with  a  military-specific 
corpus,  we  also  would  like  to  compare  the  dictionary  with  a  general  corpus,  such  as  English 
GigaWord.  This  process  would  not  validate  terms,  but  rather  would  propose  a  possible  list  for 
exclusion.  By  processing  the  dictionary  with  a  general  corpus,  we  would  be  able  to  eliminate 
general  terms,  but  also  single-word  terms  that  occur  frequently  in  both  a  general  corpus  and  a 
military-corpus.  These  unigrams  must  be  verfied  with  a  human  ground  truth  because  of  the 
appearance  of  homographs,  as  mentioned  earlier. 

The  third  proposed  method  that  we  plan  to  execute  following  this  paper  is  IDF.  The  problem 
with  TF  measurements  is  that  all  documents  and  expressions  are  considered  equally  important  in 
terms  of  assessing  relevancy.  IDF  works  to  solve  this  problem  along  with  TF  by  statistically 
identifying  how  important  a  word  is  to  a  corpus.  If  the  TF-IDF  is  high,  it  indicates  a  rare  tenn;  it 
is  considered  low  when  terms  occur  frequently. 


7.  Conclusion 


As  of  the  moment,  we  have  46.70%  of  the  dictionary  accounted  for  as  a  result  of  the  TF/term 
extraction  methods,  as  well  as  a  portion  dedicated  to  Zipf  s  Law  (8.27%).  After  all  the  previously 
mentioned  methods  have  been  executed,  we  hope  to  have  a  refined,  efficient  dictionary  that  will 
be  useful  in  the  field  as  well  as  for  more  computational  research. 
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Abstract 


When  a  solution  of  sodium  perchlorate  (NaClCfi)  in  methanol  is  drop-cast  onto  porous  silicon 
(PS)  and  dried  under  nitrogen  (N2),  an  explosive  reaction  can  be  initiated  by  heat,  friction,  or 
electrical  spark.  As  a  means  to  yield  high-performing  energetic  porous  silicon  on-chip,  a 
galvanic  corrosion  process  is  demonstrated.  This  is  a  major  advancement  for  energetic  porous 
silicon,  offering  a  more  flexible  approach  to  integration  than  previous  work  using 
electrochemical  etching.  The  new  galvanic  process  has  no  requirement  for  electrical  connection 
to  the  wafer  during  pore  formation  and  is  compatible  with  many  more  materials  than  the 
electrochemical  etch  process.  Here,  PS  films  up  to  150  pm  thick  are  fabricated  on  p-type  Si 
wafers.  Bomb  calorimetry  reveals  PS-NaCICfi  to  be  a  fuel-rich  energetic  material  with  a  heat  of 
reaction  below  the  theoretical  value.  The  reaction  velocity  is  studied  by  pattering  thin-film 
resistor  wires  across  PS  strips  3  mm  wide  and  up  to  75  mm  in  length.  An  oscilloscope  monitors 
the  wires  during  the  energetic  reaction,  and  the  time  to  failure  is  used  to  determine  the  speed  of 
the  reaction.  The  results  indicate  that  velocity  is  on  average  2130  m/s  but  can  fluctuate  hundreds 
of  m/s  along  the  length  of  the  chip. 
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1.  Introduction/Background 


The  development  and  tailoring  of  new  energetic  materials  is  critical  to  the  success  of  explosives, 
propellants,  micro-actuation,  and  power  for  next  generation  Department  of  Defense  (DoD) 
weapons  and  microsystems.  In  particular,  a  material  that  is  compatible  with  and  integrated  into 
microelectromechanical  systems  (MEMS)  could  greatly  advance  the  current  state  of  the  art  in 
MEMS  fuzing,  propulsion,  and  power  applications.  Advances  in  the  relatively  young  field  of 
nanoenergetics  have  begun  to  make  this  on-chip  integration  a  reality  (7).  Nanothermites  and 
energetic  porous  silicon  (PS)  are  the  two  main  focus  areas  in  the  field  of  nanoenergetics. 

Nano  thermite  work  investigates  the  effects  of  scaling  traditional  thermite  materials,  generally 
aluminum  and  an  oxidized  metal,  such  as  iron  oxide,  to  the  nanoscale.  Current  studies  have 
shown  improved  reaction  rate — i.e.,  increased  bum  rates — compared  to  macroscale 
thermite  materials  (2,  3).  Additionally,  a  micro-heater  has  been  realized  using  a  nanothermite 
approach  ( 4 ). 

The  work  in  this  report  focuses  on  energetic  PS,  a  widely  studied  material  that  holds  great 
potential  in  the  realization  of  novel  MEMS,  including  chemical  and  biological  sensors,  and 
optoelectronics.  When  impregnated  with  any  one  of  a  number  of  oxidizers,  PS  can  function  as 
an  energetic  material.  Because  silicon  (Si)  is  widely  used  in  micromachining  and  complimentary 
metal-oxide-semiconductor  (CMOS)  processes,  PS  can  be  integrated  on-chip  alongside  a  MEMS 
sensor.  This  new  class  of  energetic  material  can  be  used  for  applications  requiring  on-chip 
power,  propulsion,  and  fuzing.  The  strength  of  the  energetic  reaction  is  controlled  by  altering 
processing  parameters  to  tune  PS  thickness,  porosity,  specific  surface  area,  and  surface 
terminations. 

In  general,  PS  begins  as  a  single-  or  polycrystalline  Si  wafer  and  is  anodically  etched  in  a 
hydrofluoric  (HF)  acid  environment  by  applying  an  electrical  bias  to  the  Si.  The  end  result  is  a 
film  that  is  a  few  hundred  nanometers  to  several  hundred  microns  thick  and  composed  of  a 
network  of  pores  ranging  from  a  few  nanometers  to  several  microns  in  diameter.  Alternatively, 
the  film  can  be  classified  as  being  composed  of  Si  nanocrystals.  The  specific  surface  area  of  the 
material  can  be  as  high  as  1000  m2/cm\ 

For  energetic  PS,  the  pores  are  filled  with  a  strong  oxidizer,  such  as  a  perchlorate  salt  dissolved 
in  methanol  or  molten  sulfur.  After  the  methanol  evaporates  or  the  sulfur  dries,  the  pores  are 
filled  with  the  solid  oxidizer.  By  applying  heat,  mechanical  force,  or  a  spark,  the  material  readily 
ignites  with  a  bright  flash  of  light  and  a  loud  bang  approaching  120  decibels  (dB),  for  only  a  few 
milligrams  of  material. 

Early  reports  on  energetic  PS  were  primarily  qualitative  in  nature  and  proposed  reaction 
mechanisms — i.e.,  the  qualitative  strength  of  the  reaction  (as  judged  by  flame  size  and  how  loud 
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of  bang  is  produced)  with  several  solid  oxidizers  (5,  6).  One  report  proposes  an  optimal  pore 
size  for  the  most  energetic  reaction  (5).  More  recently,  attempts  at  differential  scanning 
calorimetry  (7),  traditional  calorimetry  (8),  and  Fourier  transfonn  infrared  spectroscopy  (FTIR) 
(9)  have  been  made.  Additionally,  one  report  presents  work  on  nanoporous  Si  particles  of 
several  microns  in  diameter  and  reports  ignition  by  partial  filling  of  pores  with  CuO 
nanoparticles  (10).  Most  recently,  from  the  U.S.  Army  Research  Laboratory  (ARL),  on-chip 
integration  of  energetic  PS  with  a  MEMS  sensor  is  reported  (11). 

In  general,  however,  there  is  still  very  little  quantitative  work  documenting  energetic  PS,  and  it  is 
quite  possible  that  the  full  potential  of  this  material  is  not  being  realized.  There  are  thousands  of 
journal  articles  that  document  PS  (interest  peaked  when  it  was  found  that  PS  photoluminesces), 
and  these  papers  highlight  dozens  of  characterization  techniques,  etch  recipes,  and  material 
properties.  The  nanoenergetics  field  has  to  yet  to  truly  apply  this  knowledge  base  to  energetic 
PS.  To  this  end,  a  focused  research  effort  is  being  conducted  at  ARL  to  characterize 
nanoenergetic  PS  with  a  complete  suite  of  investigations. 

With  the  knowledge  gained  from  the  nanonergetic  PS  characterization,  an  on-chip  integration  of 
the  material  can  be  realized.  Current  work  from  ARL  uses  a  process  that  can  result  in  damage  to 
the  PS  and  other  materials  on  the  wafer  during  device  fabrication.  As  a  possible  way  of  avoiding 
this  damage,  a  galvanic  cell  composed  of  Si  and  platinum  (Pt)  in  a  hydrofluoric  acid-based 
electrolyte  is  used  to  generate  PS.  While  the  galvanic  cell  method  of  PS  formation  is  not  new, 
we  have  been  able  to  create  films  up  to  50  microns  thick  that  are  stable,  highly  compatible  with 
device  fabrication  steps,  and  suggest  an  energetic  reaction  that  proceeds  faster  than  in  PS 
generated  from  a  standard  anodic  etch  process.  This  is  also  the  first  report  of  thick  (tens  of 
microns)  galvanic  PS  in  10  ohm-cm  Si  wafers.  We  duplicated  work  that  generated  thick  PS 
films  in  0.01  ohm-cm  wafers  (12),  but  these  show  non-uniform  etch  rates  across  the  exposed  Si 
surface  and  yield  qualitatively  slower  energetic  reaction  rates  compared  to  the  10  ohm-cm 
wafers. 


2.  Experiment/Calculations 


Briefly,  PS  is  typically  made  via  an  anodic  etching  process  in  a  hydrofluoric  (FIF)  acid-based 
electrolyte  (13).  A  doped  Si  wafer  of  n-  or  p-type  serves  as  the  anode,  and  a  metal  electrode  that 
is  inert  (typically  gold  (Au)  or  Pt)  in  an  HF  environment  functions  as  the  cathode.  The  Si  and 
metal  are  electrically  connected,  and  a  current  density  ranging  from  a  few  mA/cm  to  several 
hundred  mA/cm2  is  passed  such  that  oxidation  occurs  at  the  Si  surface.  As  the  etch  proceeds, 
pores  ranging  from  a  few  nanometers  to  hundreds  of  nanometers  are  produced.  The  thickness 
and  pore  size  of  the  PS  layer  is  controlled  by  current  density,  electrolyte  composition,  the  dopant 
type,  the  Si  resistivity,  and  the  etch  time.  Figure  1  shows  the  setup  of  a  typical  etch  cell,  in 
which  a  silicon  wafer  is  sandwiched  between  a  Teflon  compartment  containing  an  HF-based 
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electrolyte  and  a  bottom  block  of  Teflon.  The  o-ring  serves  to  prevent  leakage  of  the  electrolyte. 
Two  typical  optical  images  of  PS  (1  cm  area)  formed  in  an  etch  cell  are  shown.  The  conditions 
used  in  this  report  typically  result  in  a  black,  reddish  brown,  or  blue  (not  shown)  PS. 


Figure  1.  (L)  Typical  etch  cell  composed  of  Teflon  electrolyte  chamber,  a 
Au  wire  cathode,  silicon  anode,  o-ring,  and  powersupply.  Two 
optical  images  of  PS  are  included.  (R)  Galvanic  cell  arrangement 
showing  ionic  charges  (reduction  at  the  Pt  cathode,  oxidation  at  the 
Si  anode).  ShN4  acts  as  the  etch  mask  in  HF. 

Alternatively,  PS  can  be  generated  using  a  galvanic  process  that  results  in  autonomous  current 
generation  between  electrically  connected  Si  and  an  HF-resistant  metal  (Au  or  Pt),  as  shown  in 
figure  1  (14,  12).  In  this  process,  no  power  supply  is  needed.  The  open  circuit  potential  (OCP) 
of  Si  and  Pt  differ  in  an  HF-based  electrolyte.  When  electrically  connected,  Pt  raises  the  OCP  of 
Si,  and  it  has  been  shown  by  electrochemical  techniques  that  the  dissolution  of  Si  (forming  of 
pores)  is  realized  (14).  Reduction  of  the  oxidizing  agent  (either  dissolved  oxygen  in  solution  or 
addition  of  a  strong  oxidizer,  such  as  hydrogen  peroxide)  at  the  metal  cathode  drives  the  current. 
The  autonomously  generated  current  is  affected  by  electrolyte  composition  and  surface  area  ratio 
(SAR)  of  metal  to  exposed  Si. 

In  this  work,  we  use  sputtered  Pt  directly  on  a  1-10  ohm-cm  resistive  wafer  as  the  metal  contact. 
Pt  has  a  much  larger  reduction  current  than  does  Au  when  FECE  is  used  as  an  additive. 
Additionally,  we  have  found  the  Pt  can  withstand  over  30  min  of  immersion  in  HF/EtOH/FFCh. 
Ethanol  is  used  as  a  surfactant  to  aid  in  wetting  of  the  Si.  Other  surfactants,  such  as  the 
commercially  available  Triton  X-100  that  has  been  used  in  other  works  (15,  14,  16),  were  tested 
with  some  success.  Patterned  Si  wafers  using  a  Si3N4  mask  and  a  sputtered  Pt  film  on  the  back 
of  the  wafer  were  used  to  fabricate  2  mm  diameter  PS  devices.  In  the  case  of  the  2  mm  devices, 
the  SAR  is  about  2.4,  whereas  a  single  crystal  wafer  with  a  Pt  film  on  the  backside  has  a  1 .0 
SAR. 

FTIR  (Nicolet,  Thermo-Fisher  Scientific)  is  used  to  monitor  changes  in  PS  composition  as  a 
result  of  different  preparation  methods  or  treatments.  We  use  FTIR  in  reflection  mode  with  a 
diamond  window,  attenuated  total  reflection  (ATR)  fixture.  Specifically,  we  study  the  role  of 
oxidation  in  an  O2  atmosphere  at  250  °C  for  60  s  in  a  rapid  thermal  anneal  (RTA)  tool. 


67 


Calorimeter  data  (Parr  Instruments,  1109  semi-micro  bomb),  in  conjunction  with  pressure 
measurements  of  the  evolved  gas  from  the  energetic  PS  reaction,  is  presented,  revealing  how 
energy  output  and  gas  generation  change  with  Si  oxidation  prior  to  ignition. 


3.  Results  and  Discussion 


3.1  FTIR  Characterization  of  PS 

The  actual  cause  of  the  rapid  exothermic  reaction  of  PS  with  a  solid  oxidizer  is  thought  to  result 
from  a  region  of  bare  Si  nanocrystal  becoming  exposed  to  the  oxidizer.  The  high  concentration 
of  a  strong  oxidizer  in  the  presence  of  exposed  Si,  which  has  a  higher  exothermic  energy  release 
after  oxidation  than  carbon,  leads  to  a  rapid  chain  reaction  with  a  large  energy  release  occurring 
in  milliseconds.  “Fresh”  PS — PS  only  several  hours  old  or  less — exhibits  a  hydrogen-terminated 
Si  surface  that  prevents  the  slow  oxidation  of  the  Si  nanocrystals  in  the  ambient  environment.  It 
is  also  thought  that  the  hydrogen  is  what  leads  to  the  gas  evolution  during  ignition  of  the 
nanoenergetic  PS.  From  rudimentary  tests  in  which  a  mass  is  propelled  from  the  surface  of  a  PS 
ignition,  it  is  evident  that  gas  is  being  generated  from  the  reaction.  Equation  1  is  proposed  as  the 
reaction  equation  for  nanoenergetic  PS  without  consideration  of  any  other  species,  such  as 
hydrogen  or  crystalline  water,  in  the  oxidizer.  Since  the  reaction  products  of  the  fuel,  Si,  and  the 
oxidizer,  NaClO/t,  are  all  solid,  the  gas  evolution  must  be  coming  from  another  source. 

2Si  +  NaC104  -»  2Si02  +  NaCl  (1) 

Generation  of  hydrogen  gas  and  water  vapor  as  a  result  of  the  hydrogen  surface  termination 
being  driven  is  likely  (5).  One  source  suggests  that  a  mild  oxidation  of  PS  does  not  diminish  the 
exothermic  energy  release  after  ignition;  it  may  improve  stability  of  the  removing  relatively 
weak  Si-H  bonds  by  either  forming  O-Si-H  or  removing  all  of  the  Si-H  to  form  a  thin  Si02  layer 
(US  patent:  6984274).  Figure  2  shows  FTIR  data  of  “fresh”  and  oxidized  (details  described  in 
the  experimental  section)  Si.  Two  variants  of  PS  Si  are  shown:  (1)  PS  from  an  electrochemical 
etch  at  20  mA/cm2  for  30  min  in  1 : 1  (49%  HF:EtOH  by  volume),  and  (2)  PS  from  a  galvanic 
process  etched  for  20  min  in  1:11 :22  (70%  H202:Et0H:49%  HF  by  volume),  with  a  1 .0  SAR  of 
Pt  to  Si. 

The  relative  peak  intensities  in  figure  2  reveal  that  the  oxidation  performed  in  relatively  mild 
conditions  and  short  time  does  not  completely  remove  the  hydrogen  termination,  but  there  is  a 
clear  increase  in  Si-0  bonds,  indicative  of  a  partial  oxidation  of  the  surface.  It  is  also  noted  that 
the  susceptibility  of  the  PS  to  accidental  ignition  by  mechanical  stimulation  did  decrease  as 
suggested  in  (US  patent:  6984274).  The  two  methods  of  PS  preparation  yield  similar  FTIR  data, 
with  the  galvanic  process  showing  a  slight  resistance  to  oxidation,  as  the  height  of  the  peaks  at 
2100  cnT1  did  not  diminish  as  much  as  the  electrochemical  etch. 
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Figure  2.  FTIR  of  “fresh”  and  oxidized  PS  from  a  standard  electrochemical  etch  and  a 
galvanic  etch  process. 

3.2  Calorimetry  of  PS 

To  determine  the  effect  of  oxidation  on  energetic  PS,  calorimetric  measurements  of  “fresh”  and 
oxidized  PS  with  NaClC>4  oxidizer  were  made.  In  each  case,  the  pressure  evolved  from  each 
reaction  was  measured  in  a  separate  experiment  using  a  standard  pressure  gauge  connected  to  the 
Parr  calorimetry  bomb.  Figure  3  shows  the  temperature  rise  as  a  function  of  time;  both  variants 
of  PS  show  similar  behavior.  The  exothermic  energy  from  each  reaction  is  very  similar  at  5.4 
kJ/g  for  oxidized  PS  and  5.7  kJ/g  for  “fresh”  PS. 


Figure  3.  Temperature  rise  during  bomb  calorimetery. 

The  two  variants  differ  dramatically  in  gas  volume  generated.  “Fresh”  PS  yields  about 
0.0129mol/g,  whereas  the  oxidized  Si  generates  about  0.004  mol/g.  This  suggests  that  the 
hydrogen  termination  plays  a  significant  role  in  gas  generation,  but  minimally  affects  the 
exothermic  energy  of  the  reaction.  This  discovery  is  likely  to  provide  avenues  for  “tuning”  the 
gas  generation  of  the  reaction  for  different  applications.  Further  analysis  and  testing  of  additional 
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PS  variants  are  required  to  fully  understand  this  result,  but  initial  propagation  rate  studies,  from 
high  speed  video  analysis,  indicate  that  the  oxidized  sample  reacts  slower  than  the  “fresh” 
sample. 

Additionally,  higher  propagation  speeds  can  be  correlated  to  louder  energetic  reactions  (120 
db+),  and  both  the  “fresh”  and  galvanic  PS  achieve  the  loudest  ignitions.  The  galvanic  PS 
especially  is  highly  reproducible  in  ignition  sound,  with  even  the  2  mm  devices  producing  a  large 
flash  and  loud  bang. 


3.3  On-chip,  Nanoenergetic  Galvanic  PS 


Energetic  PS  devices  are  currently  realized  by  using  a  SfiN4  etch  mask  to  pattern  a  Si  wafer. 
Next,  an  electrochemical  etch  fabricates  PS  in  the  exposed  regions  of  the  Si3N4  mask.  Finally,  a 
metal  bridgewire  (shown  in  figure  4a)  is  patterned  on  the  PS.  A  current  can  be  passed  through 
the  bridgewire  to  cause  localized  heating  of  the  PS  and  subsequent  ignition  when  oxidizer  is 
present.  During  the  processing  steps  to  realize  these  devices,  the  PS  is  degraded  and  large 
cracks  are  formed  that  can  result  in  bridgewire  damage  (figure  4b).  Additionally,  the  etch  depths 
are  limited  by  this  process. 


Figure  4.  SEM  images,  (a)  Bridgewire  deposited  on  bare  Si;  inset)  expanded 
view  of  bridgewire.  (b)  Bridgewire  deposited  on  electrochemically 
etched  PS.  (c,d)  Bridgewire  after  galvanic  PS  generation.  The  higher 
H2O2  concentration  in  (c)  results  in  a  more  aggressive  etch  than  an  (d). 

In  an  effort  to  realize  structurally  stable  PS  and  more  reliable  bridgewire  structures,  as  well  as  to 
reduce  the  number  of  processing  steps  and  increase  device  yield  per  wafer,  galvanic  PS 
generation  was  used  to  create  PS  films  on  the  order  of  20-50  pm  thick  with  limited  bridgewire 
degradation.  In  this  process,  the  bridgewire  is  patterned  on  bare  Si  that  is  then  galvanically 
etched  to  generate  PS.  Results  indicate  that  the  bridgewire  remains  intact,  and  the  PS  is 
structurally  stable  with  limited  cracking  present. 
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Figure  5  presents  SEM  cross  section  images  of  galvanic  PS  with  different  SAR  and  different 
electrolyte  compositions.  The  etch  depth  increases  both  with  increasing  H2O2  concentration,  as 
well  as  SAR.  Additionally,  HF  concentration  and  EtOH  concentration  affect  the  structural 
integrity  of  the  film.  The  wafer  resistivity  also  drastically  affects  the  PS  film  quality  and  etch 
depth.  More  results  are  needed  to  develop  a  more  exact  correlation  of  SAR  and  electrolyte 
composition  to  film  thickness  and  quality,  but  these  preliminary  results  indicate  that  it  is  possible 
to  generate  stable  porous  Si  layers  at  least  50  pm  thick  using  the  galvanic  process. 


Figure  5.  SEM  cross  section  images  of  galvanic  PS  with  different 
SAR  and  different  electrolyte  compositions. 

Currently,  appropriate  galvanic  etch  conditions  have  been  used  to  fabricate  stable  PS  films  after 
bridgewire  deposition.  These  devices  have  been  wirebonded  such  that  ignition  can  be  tested. 
While  no  quantitative  data  has  been  attained,  the  reaction  propagation  velocity  in  small  devices 
(2  mm  diameter)  appears  to  be  much  greater  than  the  standard  electrochemical  PS  devices. 
Additionally,  the  bridgewire  is  of  high  quality  and  less  prone  to  failure  due  to  cracking  of  the 
porous  film.  This  process  is  still  very  young  but  could  realize  very  unique  device  architectures, 
faster  device  fabrication,  and  more  reproducible  ignitions. 

3.4  On-chip  Velocity  Measurements 

The  velocity  at  which  the  energetic  reaction  between  PS  and  NaC104  proceeds  is  an  important 
parameter.  The  velocity  indicates  if  the  material  is  truly  detonating  or  is  slowly  burning. 
Molecular  energetic  materials  contain  the  fuel  and  oxidizer  in  one  molecule  and,  thus,  oxidation 
occurs  extremely  quickly.  Reaction  velocities  of  these  materials  approach  10  km/s.  The  PS- 
NaC104  is  not  expected  to  reach  this  velocity  since  the  reaction  involves  the  decomposition  of 
NaC104  prior  to  oxidizing  PS.  However,  modifying  the  PS  morphology,  either  through  etching 
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conditions  or  with  heat  and  chemical  treatments,  is  expected  to  change  the  velocity  of  the 
reaction. 


High  velocity,  explosive  reactions  are  often  studied  with  a  high  speed  camera,  but  this  is  a  time¬ 
intensive  and  costly  technique.  For  an  on-chip  measurement,  resistive  wires  that  are  very  similar 
to  the  initiator  wire  previously  discussed  are  lithographically  patterned  on  the  PS.  The  reaction 
is  initiated  at  one  of  the  chip,  and  as  it  progresses,  it  breaks  the  wires,  as  shown  schematically  in 
figure  6.  Initially,  9  V  is  applied  across  the  wires,  and  this  voltage  is  monitored  with  an 
oscilloscope.  The  voltage  drop  at  each  wire  is  broken. 

9V  Reaction  direction 


Initiateon 

thisend 


Figure  6.  Schematic  of  on-chip  velocity  measurement. 

A  representative  plot  of  the  wires  failing  is  shown  in  figure  7.  The  reason  for  the  increase  in  the 
voltage  after  the  wires  fail  before  once  again  settling  to  0  V  could  relate  to  either  properties  of 
the  oscilloscope  or  effects  from  the  energetic  reaction  reconnecting  the  wires  via  conductive 
plasma.  For  a  SAR  of  5. 2-6.4,  the  velocity  is,  on  average,  2130  m/s.  This  velocity  can  fluctuate 
hundreds  of  m/s  along  the  chip,  however.  This  may  relate  to  areas  where  NaCICfi  dries 
differently  or  is  not  reacting  with  the  PS  in  the  same  manner  as  other  regions.  Initial  work 
reveals  that  when  low  concentration  NaCICfi  is  used  (0.8M-1.6M),  the  reaction  velocity  slows  to 
hundreds  of  m/s  and,  in  some  cases,  tens  of  m/s.  Additionally,  high  speed  video  analysis  should 
be  conducted  alongside  these  measurements  to  verify  some  of  the  findings,  although  the  average 
velocity  presented  here  is  within  an  expected  range  (9).  Alongside  calorimetric  data,  this 
velocity  test  can  be  performed  to  gain  insight  into  how  PS  morphology,  surface  termination,  or 
use  of  an  oxidizer  other  than  NaClCfi  affects  the  reaction. 
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Figure  7.  Representative  velocity  data  from  the  on-chip 
measurement  with  wire  spacing  of  5  mm. 


4.  Summary  and  Conclusions 


Nanoenergetics  is  a  very  young  field  of  research,  and  limited  quantitative  data  and  material 
characterization  is  presented  in  the  literature.  The  focus  of  this  report  is  nanoenergetic  PS. 
Calorimetry,  FTIR,  SEM,  and  velocity  data  are  presented  for  several  variants  of  PS.  Quantitative 
data  will  aid  in  the  ability  to  fine-tune  the  PS  morphology  to  yield  desired  gas  volume, 
propagation  speed,  and  energy  density.  A  new  galvanic  process  to  form  energetic  PS  is 
introduced,  and  initially  suggests  an  easier,  faster,  and  more  reliable  route  to  on-chip  integration. 
Additional  benefits  seem  to  include  increased  propagation  rate  and  more  reliable  bridgewire 
structures.  Future  work  will  include  further  study  into  the  structure  of  PS  and  how  it  affects 
velocity  and  energy  yield. 
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Abstract 


The  protective  antigen  (PA)  protein  secreted  by  the  Bacillus  anthracis  bacterium  combines  with 
either  edema  factor  or  lethal  factor  protein  to  form  either  edema  toxin  or  lethal  toxin.  Since  PA  is 
necessary  to  the  formation  of  both  toxins,  it  can  be  used  as  a  target  to  detect  the  presence  of  B. 
anthracis.  Because  antibody-antigen  detection  methods  are  time  consuming  and  require  specific 
conditions  for  testing,  alternative  methods,  such  as  peptide-antigen  detection,  are  being 
considered.  Enzyme-linked  immunosorbent  assay  (ELISA)  is  commonly  used  to  test  for  binding 
affinity  of  proteins  and  enzymes,  and  employs  the  use  of  visible  detection  of  substrate-enzyme 
reactions.  Binding  affinity  is  quantified  based  on  the  dissociation  constant  (Kd)  of  a  peptide- 
antigen  or  antibody-antigen  complex,  with  low  Kd  values  representing  high  affinity.  A 
fluorogenic  substrate  was  reacted  against  a  15 -amino  acid  peptide  bound  to  a  PA-enzyme 
conjugate  using  the  ELISA  method.  The  determined  Kd  for  the  peptide  was  in  the  range  of 
680-728  nM,  while  the  binding  antibody  specific  to  PA  was  in  the  range  of  ~3— 5  nM.  Truncated 
segments  of  the  full-length  peptide  will  be  tested  to  determine  a  sequence  with  the  highest 
binding  affinity  to  PA,  providing  a  more  efficient  alternative  for  field  use. 
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1.  Introduction/Background 


Following  its  use  in  2001,  Bacillus  anthracis  and  its  associated  toxins  have  received  much 
attention.  B.  anthracis  secretes  three  distinct  proteins,  identified  collectively  as  anthrax  toxin: 
edema  factor  (EF,  89  kDa),  lethal  factor  (LF,  90  kDa),  and  protective  antigen  (PA,  83  kDa). 
Cleaving  of  PAg3  on  the  cell  surface  into  PA20  and  the  activated  PA63  protein  by  furin  proteases 
exposes  binding  sites  for  EF  and  LF  ( 1 ,  2).  PA  is  easily  passed  through  the  cellular  membrane, 
and  thus  any  bound  EF  or  LF  protein  can  be  introduced  inside  the  cell  and  undergo  its  respective 
biological  toxins — edema  toxin  (EdTx)  and  lethal  toxin  (LeTx)  (figure  1)  (2-4). 


Figure  1.  Mechanism  of  EF  and  LF  introduction  inside  cell  membrane  via  binding  with 
PA  ( 3 ). 

Biochemical  and  computational  research  addressing  different  mechanisms  of  molecular 
recognition  is  currently  being  conducted  to  gain  a  better  understanding  of  binding  affinities  and 
target  specificity  using  synthetic  alternatives.  Peptide  recognition  elements  have  emerged  as  a 
potential  synthetic  alternative  to  antibodies.  Antibody  production  against  a  specific  protein  often 
requires  in  vivo  conditions  and  several  months  of  isolation,  whereas  target-specific  peptide 
binders  can  be  sequenced  and  isolated  in  vitro  within  2-4  days  (5).  The  overall  focus  of  our 
research  is  to  develop  a  high-affinity  peptide  with  the  sequence  necessary  for  PA  binding  that 
would  provide  low-cost,  high-stability,  rapid,  and  accurate  detection  methods  for  such  virulent 
toxins  as  EdTx  and  LeTx. 


81 


My  current  project  to  this  end  goal  centers  around  the  use  of  enzyme-linked  immunosorbent 
assay  (ELISA)  to  determine  the  binding  affinity  of  the  15  amino  acid  SM545  peptide,  expressed 
in  terms  of  the  dissociation  constant,  Kd  (units  of  nM).  For  the  equilibrium  dissociation  reaction 
of  a  PA-peptide  complex, 

PA*  Peptide  <->  PA  +  Peptide.  (1) 

The  Kd  is  determined  via  the  following  expression: 

Kd  =  [PA]  [Peptide]  /  [PA*  Peptide  complex].  (2) 

This  paper  focuses  on  ELISA  and  its  biochemical  uses  in  peptide-protein  interaction  and 
detection.  Direct  heterogeneous  ELISA  is  a  simple,  rapid,  and  sensitive  method  that  exploits  the 
target-specific  binding  behavior  of  proteins.  It  employs  the  use  of  96-well  polystyrene  plates 
(figure  1)  in  which  reagents  are  sequentially  added  and  directed  at  a  solid  phase,  and  several 
washing  steps  are  performed  to  separate  bound  (reacted)  and  unbound  (unreacted)  reagents  (6). 
Detection  is  based  on  the  chromogenic  or  chemifluorescent  interaction  between  an  enzyme- 
linked  protein  and  a  reacting  substrate  that  produces  a  readable  signal.  Because  biological 
macromolecules,  such  as  proteins,  have  the  potential  to  denature  and  lose  proper  binding 
function,  optimization  of  a  consistent  method  is  a  key  factor  in  obtaining  reproducible  results. 
Most  of  my  current  research  has  dealt  with  identifying  these  experimental  complications  and 
making  the  necessary  modifications.  Following  successful  KD  determination  of  the  full-length 
SM545  peptide,  we  will  analyze  fragmented  truncates  of  varying  lengths  to  determine  the  exact 
sequence  of  the  peptide  necessary  for  PA  binding. 


2.  Methods  and  Materials 


2.1  Procedure  for  Direct  Heterogeneous  ELISA 

The  general  procedure  of  the  direct  heterogeneous  ELISA  is  outlined  in  figure  2  (6).  The  steps 
are  as  follows: 

1 .  A  1 00- pL  sample  of  peptide  diluted  in  buffer  is  added  to  the  solid  phase,  incubated  for  2 
h  at  room  temperature  (i),  and  allowed  to  passively  adsorb  to  the  well  walls  (ii)  (6-8). 

2.  After  incubation,  the  wells  are  emptied  of  the  sample,  and  300  pL  of  a  blocking  agent 
(neutral  buffer  +  detergent)  is  added  to  the  wells  and  incubated  for  3  h.  This  blocking  step 
is  essential  to  coat  all  available  binding  sites  on  the  solid  phase  to  deter  further  passive 
adsorption  by  reaction  reagents  (6,  9,  10). 
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Figure  2.  Procedural  outlines  for  direct  heterogeneous  ELISA. 

3.  The  wells  are  emptied  again  and  flooded  three  times  with  300-pL  aliquots  of  the  blocking 
agent.  These  washing  steps  remove  any  unbound  (unreacted)  peptide. 

4.  The  antigen-labeled  enzyme  (100  pL)  is  then  added  to  the  reaction  wells  and  incubated  at 
room  temperature  for  1  h  (iii). 

5.  Two  additional  washing  steps  are  performed  using  a  neutral  buffer  to  remove  residual 
detergent  present  in  the  wells  from  the  blocking  agent  ( 6 , 11).  At  this  time,  peptides  are 
bound  to  the  solid  phase  as  well  as  to  the  enzyme-labeled  antigen.  Chromophoric  or 
fluorophoric  substrate  (100  pL)  (see  section  2.2)  is  added  to  the  wells  and  incubated  for 
~  10-30  min  (iv).  A  readable  signal  is  observed  at  this  stage.  Enzymatic  activity  is 
stopped,  and  the  colorimetric  or  chemifluorescent  reaction  discontinued  via  denaturation 
of  enzymes  ( 6 ).  The  reaction  is  quantified  using  a  MicroWell  plate  reader  (vii). 
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2.2  Reaction  Reagents 


The  purpose  of  the  following  experiment  was  to  investigate  the  binding  affinity  of  the 
polystyrene-tagged  SM545  peptide  (molecular  weight  [MW]:  3515  g/mol)  dissolved  in  a 
degassed  0.2-M  carbonate/bicarbonate  buffer  with  a  pH  =  9.4.  PA  mAb  (MW:  150,000  g/mol) 
antibody  dissolved  in  the  same  buffer  was  used  as  a  positive  control.  The  carbonate/bicarbonate 
buffer  was  chosen  as  a  solvent  because  it  contained  no  competing  proteins  and  also  served  as  a 
negative  control  to  ensure  that  unreacted  reagents  were  not  present  in  the  final  quantified  sample. 
Phosphate-buffered  saline  (PBS)  +0.1%  Tween-20  was  used  to  block  and  wash  the  wells.  The 
0.5  ug/ml  PA  conjugated  with  horseradish  peroxidase  (HRP)  was  diluted  in  PBS  and  used  as  the 
enzyme-linked  antigen.  HRP  reacts  via  oxidation  with  chromogenic  substrates,  such  as  3,  3’,  5, 
5’-  tetramentylbenzidine  (TMB),  or  chemifluorescent  substrates,  such  as  10-acetyl-3-7- 
dihydroxyphenoxazine  (ADHP).  A  blue  product  is  observed  when  the  HRP-conjugated  antigen 
reacts  with  TMB,  while  a  resorufin,  a  fluorescent  product,  is  observed  when  reacted  with 
ADHP.  The  TMB  reaction  is  stopped  using  a  high  concentration  of  a  strong  acid  such  as  0.5-M 
H2SO4  or  1-M  H2SO4  in  8-M  acetic  acid,  yielding  a  yellow  end  product  with  absorbance  at 
450  nm.  The  HRP -ADHP  complex  fluoresces  at  570  nm. 

2.3  Plate  Setup  and  Sample  Prep 

The  Nunc  Immuno  96  Microwell  plates  were  used  for  the  solid  phase.  The  plate  was  set  up  as 
shown  in  figure  3.  The  edge  plates  were  left  without  samples  to  avoid  temperature  differences 
and  edge  effects  (see  section  2.4)  (6,  12).  The  concentrations  of  the  samples  in  the  wells  are 
given  in  table  1 . 


Figure  3.  Pictorial  representation  of  a  96-well  plate  setup. 
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Table  1.  Antibody  and  peptide  sample  concentrations  (ug/mL). 


5.0  3.0  2.5  2.25  2.0  1.75  1.5 

pg/ml 

1.25  1.0  0.50 

5.0  3.0  2.5  2.25  2.0  1.75  1.5 

pg/ml 

1.25  1.0  0.50 

5.0  3.0  2.5  2.25  2.0  1.75  1.5 

pg/ml 

1.25  1.0  0.50 

5.0  3.0  2.5  2.25  2.0  1.75  1.5 

pg/ml 

1.25  1.0  0.50 

2.0  1.75  1.5  1.25  1.0  0.75  0.50 

pg/ml 

0.25  0.125  0.06 

0.2-M  Carbonate/Bicarbonate  Buffer 

2.4  ELISA  Complications  and  Modifications 

In  the  ELISA  method,  the  use  of  peptides  and  protein  makes  several  factors  important  to 
successful  and  repeatable  data  collection  and  Kd  determination.  The  first  factor  is  amount.  The 
same  aliquot  of  liquid  reagent  needs  to  be  delivered  to  all  wells  in  each  step  of  the  procedure.  To 
offset  this  factor,  the  same  calibrated  pipettes  were  used  in  all  trials  to  maintain  instrumental 
precision,  and  the  delivery  method  was  standardized  across  all  wells  in  the  plate  to  minimize 
human  error.  During  all  incubation  steps,  the  reaction  plates  were  placed  on  a  rotator  to  ensure 
even  distribution  of  reagents  within  the  wells.  With  the  constant  and  sequential  emptying  and 
flooding  of  the  high-proximity  wells  with  liquids,  the  possibility  of  cross  contamination  arises. 
Precautions,  such  as  constant  changing  of  pipette  tips  and  removal  of  liquid  via  pipette  withdraw 
instead  of  a  simple  overturn  of  reaction  plate,  were  taken  to  ensure  contaminations  were 
minimized.  Lastly,  optimization  of  molecules,  such  as  peptides  and  proteins  (antigens  and 
antibodies),  is  important  to  maintain  working  reactants.  To  counteract  the  potential  of  light  and 
temperature-induced  denaturation,  aluminum  foil  was  placed  over  the  reaction  plate  during  all 
incubation  periods,  and  all  peptide  and  protein  samples  were  kept  cold  when  not  in  use. 
Temperature  differences  across  the  well  are  also  inherently  present  in  the  plates  themselves. 
Peripheral  wells  have  a  higher  evaporation  rate  than  the  interior  well  and  therefore  are  cooler  (7). 
Lor  this  reason  the  outer  wells  of  the  plate  were  not  used,  and  evaporation  rates  were  slowed 
down  by  sealing  the  plate  and  lid  cover  together  with  Parafilm  during  incubation  steps. 
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3.  Results 


Several  trials  of  the  full-length  SM545  peptide  and  the  PA  mAh  antibody  ELISA  procedure  were 
performed,  and  the  degree  of  dissociation  (KD)  was  determined  by  plotting  normalized 
fluorescence  vs.  molar  concentration  of  peptide  and  antibody  and  fitted  to  a  sigmoid  function 
(figure  4).  The  molecular  weight  of  the  full-length  polystyrene-tagged  peptide  is  35 15  g/mol  and 
that  of  the  antibody  is  150,000  g/mol.  From  this,  the  following  molar  concentrations  of  peptide 
and  antibody  can  be  calculated  (table  2). 


Figure  4.  KD  determination  curves  of  SM545-15  peptide  and  mAb  using  ADHP. 
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Peptide  ( rows  B&  C):  MW=  3515g/mol 

Well 

number 

Volume-based 

concentration 

Molar 

concentration 

2 

5.00  pg/ml 

1 .42  x10"b  M 

3 

3.00  pg/ml 

8.53  x10_/  M 

4 

2.50  pg/ml 

7.11  x10_/  M 

5 

2.25  pg/ml 

6.40  x10_/  M 

6 

2.00  pg/ml 

5.69  x10"'  M 

7 

1.75  pg/ml 

4.98  x10_/  M 

8 

1.50  pg/ml 

4.27  x10-'  M 

9 

1.25  pg/ml 

3.56  x10_/  M 

10 

1.00  pg/ml 

2.84  x10-'  M 

11 

0.50  pg/ml 

1 .42  x10_/  M 

Antibody  (  Row  E)  MW=  150,000g/mol 

Well 

Volume-based 

Molar 

number 

concentration 

concentration 

2 

2.00  pg/ml 

1.33  xl  0-8  M 

3 

1.75  pg/ml 

1.17x1 0-8  M 

4 

1.50  pg/ml 

I.OOx  10-8M 

5 

1.25  pg/ml 

8.33x  10-9M 

6 

1.00  pg/ml 

6.67x  10-9M 

7 

0.75  pg/ml 

5.00x  10-9M 

8 

0.50  pg/ml 

3.33  x10-9M 

9 

0.25  pg/ml 

1.67x  10-9  M 

10 

0.125  pg/ml 

8.33x10-10  M 

11 

0.060  pg/ml 

4.00x10-10  M 

Table  2.  Converted  molar  concentrations  of  antibody  and  peptide  samples. 

From  this  data,  we  find  that  the  Kd  of  the  full-length  SM545  peptide  is  704.3  ±  24.7  nM  ±  while 
the  mAh  antibody  has  an  equilibrium  dissociation  constant  of  4.90  ±  0. 14  nM.  A  follow-up  to  the 
full-length  SM545  peptide  was  to  analyze,  using  fluorescence,  a  5-amino  acid  truncated  segment 
of  the  SM545  peptide  with  molecular  weight  of  2338  g/mol  (figure  5).  The  same  volume-based 
peptide  concentrations  as  presented  in  table  2  were  used. 


Figure  5.  KD  determination  curves  of  SM545-5  truncated  peptide  and  mAb  antibody  using  ADHP. 
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After  eight  trials,  the  Kd  of  the  truncated  peptide  was  determined  to  be  467.6  ±  178.0  nM.  As 
shown  in  figure  4,  the  data  does  not  fully  fit  a  sigmoidal  function.  The  lower  limit  of  the  peptide 
concentration  will  be  extended  in  order  to  provide  a  more  accurate  determination  of  the  truncated 
peptide’s  dissociation  constant. 


4.  Summary  and  Conclusions 


From  these  results,  the  full-length  SM545  peptide  binds  with  lower  affinity  to  PA  than  the  mAh. 
Preliminary  results  also  suggest  that  although  the  SM545-5  truncate  also  binds  with  lower 
affinity  to  PA  than  mAh,  it  binds  at  a  higher  affinity  than  the  full-length  peptide.  For  this  reason, 
experiments  will  continue  using  truncated  peptide  fragments  of  varying  lengths  of  the  peptide  to 
determine  the  sequence  with  the  highest  binding  affinity  to  PA.  Although  it  is  feasible  to  use 
TMB  for  signal  detection,  analysis  of  fluorescence  of  the  truncated  peptides  using  ADHP  will  be 
used  exclusively  as  the  detection  method  for  remaining  experimentation  because  it  provides 
easily  reproducible  results. 
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Abstract 


The  analysis  of  fragments  removed  from  Soldiers  killed  in  action  (KIA)  is  vital  to  understanding 
the  vulnerabilities  to  threats  that  they  face  in-theater,  as  well  as  for  understanding  enemy  tactics, 
techniques,  and  procedures  (TTPs).  As  part  of  the  Joint  Trauma  Analysis  and  Prevention  of 
Injury  in  Combat  (JTAPIC)  program,  fragments  of  particular  interest  are  removed  during 
autopsy  by  the  Office  of  the  Armed  Forces  Medical  Examiner  (OAFME)  and  sent  to  the  U.S. 
Army  Research  Laboratory’s  (ARL)  Warfighter  Survivability  Branch  (WSB)  for  analysis.  The 
fragments  are  processed  and  analyzed,  and  these  results  are  entered  into  the  JTAPIC  Fragment 
and  Material  Database.  The  processing  and  analysis  data  includes  fragment  photographs, 
dimensions,  weight,  density,  three-dimensional  scans,  and  elemental  analysis  results.  The 
JTAPIC  fragment  analyses  enable  analysts  to  determine  the  identification  of  the  fragment  and 
determine  the  fragment's  origin.  This  information  is  used  in  event  recreations,  modeling,  and 
simulation.  This  paper  focuses  on  all  phases  of  the  fragment  analysis,  including  receipt  of 
fragments  and  coordination  with  the  Weapons  and  Materials  Research  Directorate  (WMRD)  for 
scanning  and  sterilization.  Upon  completion  of  WMRD’s  efforts,  the  students  processed  and 
analyzed  all  fragments  and  entered  their  results  into  the  database. 
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1.  Introduction/Background 


The  Joint  Trauma  Analysis  and  Prevention  of  Injury  in  Combat  (JTAPIC)  program  was 
established  at  the  U.S.  Army  Medical  Research  and  Material  Command  (USAMRMC),  based  on 
DODD  6025. 25E  (7).  JTAPIC  is  a  partnership  whose  purpose  is  to  collect,  integrate,  and 
analyze  operational  and  injury  data.  Its  goal  is  to  understand  the  vulnerabilities  to  threats  that 
U.S.  Soldiers  face  in-theater  as  well  as  enemy  tactics,  techniques,  and  procedures  (TTPs). 
Fragment  analysis  is  one  project  that  was  established  as  a  part  of  the  JTAPIC  Program.  During 
an  event  in-theater,  fragments  are  often  embedded  in  Soldiers.  Fragments  are  removed  during 
autopsy  by  the  Office  of  the  Armed  Forces  Medical  Examiner  (OAFME).  Analysts  within  the 
Warfighter  Survivability  Brach  (WSB)  of  the  U.S.  Army  Research  Laboratory  (ARL)  receive, 
process,  and  analyze  fragments  of  interest.  Fragments  are  analyzed  to  determine  their  elemental 
composition  and  origin.  Knowing  the  composition  of  the  fragments  can  assist  with  event 
recreation,  modeling,  and  simulation,  as  well  as  provide  insight  into  fragments  still  embedded  in 
those  wounded  in  the  same  event.  All  fragment  information  is  stored  in  the  JTAPIC  Fragment 
and  Material  Database. 

The  JTAPIC  Fragment  and  Material  Database  was  created  by  WSB  to  collect  and  store  the 
fragment  information.  The  database  provides  an  easily  accessible  storage  method  and  search 
tool  that  can  be  used  by  all  WSB  analysts  working  on  fragment-related  projects.  The  database 
contains  information  including  fragment  photographs,  dimensions,  weight,  density,  three- 
dimensional  (3-D)  scans,  and  elemental  analysis  results. 

The  fragment  processing  and  analysis  procedures  are  described.  Examples  of  the  analysis  results 
are  included  in  this  report.  A  total  of  47  fragments  from  8  cases  were  analyzed.  These  results 
were  briefed  to  the  Dismounted  Soldier  Incident  Analysis  Working  Group  and  were  provided  to 
all  JTAPIC  partners. 


2.  Experiment/Calculations 


JTAPIC  fragment  analysis  is  an  ongoing  project  run  by  WSB  analysts.  OAFME  transfers  cases 
of  interest  to  WSB  for  processing.  A  forensic  investigator  from  OAFME  transfers  custody  of 
each  fragment  delivery  to  a  WSB  evidence  custodian.  The  fragments  are  then  transported  to 
Aberdeen  Proving  Ground  (APG),  MD.  The  WSB  evidence  custodian  contacts  all  personnel 
involved  in  fragment  analysis  to  develop  an  analysis  schedule  for  each  delivery.  The  following 
organizations  are  involved  in  the  fragment  analysis  process:  the  Safety  and  Health  Physics 
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Branch  of  the  Associate  Director  for  the  Laboratory  Operations  Directorate  (ALDO)  and  the 
Multi-Functional  Materials  Branch  (MMB)  of  the  ARL  Weapons  and  Materials  Research 
Directorate  (WMRD). 

Before  WSB  can  handle  the  fragments  received  from  OAFME,  they  must  first  be  scanned  for 
radiation,  in  particular,  for  the  presence  of  depleted  uranium.  The  fragments  also  must  be 
sterilized  to  remove  any  potential  biological  hazard.  The  radiation  scanning  and  sterilization  is 
conducted  at  WMRD.  The  radiation  scanning  is  conducting  using  a  hand-held  Geiger  counter. 
Once  the  fragments  have  been  confirmed  to  be  free  of  radiation,  photographs  are  taken  of  the 
fragments  from  each  case.  The  chain  of  custody  (COC)  paperwork,  DA  Form  4173,  is  included 
in  the  photograph,  so  that  the  unique  medical  examiner  (ME)  identification  number  is  visible  as  a 
reference  within  the  photos.  These  photographs  are  vital  as  a  precautionary  measure  to  record 
the  condition  in  which  the  fragments  arrived  at  WMRD  in  case  any  of  the  fragments  should  be 
harmed  during  the  sterilization  process.  The  Fragment  Evidence  Logbook  is  used  to  record  the 
date  and  time  of  evidence  transfers.  Notes  on  any  distinctive  features  on  the  fragments  and  any 
ME  information  written  on  the  fragment  containers  (such  as  recovery  location)  is  also  recorded 
in  the  logbook  (2). 

Once  all  information  has  been  recorded  in  the  logbook,  the  fragments  can  be  sterilized.  There 
are  two  separate  sterilization  procedures  used  by  WMRD.  Metallic  fragments  are  sterilized 
using  a  VSR  Model  AS  12  autoclave.  The  autoclave  heats  the  fragments  at  100  °C  for  100  min  at 
17  psi  of  pressure.  Fragments  that  appear  to  be  plastic  are  separated  from  the  rest  of  the 
fragments  during  the  photographing  stage.  These  fragments  are  soaked  in  formalin  for  24  h, 
which  provides  the  same  level  of  sterilization  as  the  autoclave  method,  but  ensures  that  the 
plastic  will  not  melt  (2). 

After  sterilization  is  completed,  the  fragments  are  transported  by  the  WSB  evidence  custodian 
for  analysis.  The  fragments  in  each  case  are  first  placed  into  separate  containers  and  the  number 
of  fragments  for  each  case  is  noted  in  the  database  and  logbook.  Each  fragment  is  examined 
individually.  The  weight  of  each  fragment  is  taken  using  a  Scout™  Pro  SP402  scale,  which  has 
a  maximum  weight  measurement  of  400  g  with  a  readability  of  +0.01  g.  The  length,  width,  and 
depth  of  each  fragment  are  then  measured  using  a  0-200  mm  digital  caliper.  Observations  of 
any  distinguishing  coloring,  marks,  or  features  of  each  fragment  are  made  using  a  light 
microscope  and  are  recorded  in  the  database  along  with  the  weight  and  dimensions  previously 
noted.  The  fragments  are  photographed  again  from  all  sides  using  an  L-shaped  forensic  evidence 
ruler  and  grid  paper  to  indicate  the  size  of  the  fragment  in  the  photograph.  The  ME  number  is 
written  underneath  the  grid  setup  for  reference.  These  photographs  are  cropped  for  uniformity 
and  entered  into  the  database  under  the  ME  identification  number  that  corresponds  to  each 
specific  fragment.  These  photographs  are  important  because  they  preserve  the  image  of  the 
fragment,  which  could  be  used  in  a  military  court  of  law  as  evidence  of  a  homicide  (2). 
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Three-dimensional  scans  of  the  fragments  are  important  for  the  same  reason.  Any  fragment  that 
has  been  removed  directly  from  the  body  of  a  Soldier  killed  in  action  (KIA)  is  scanned  using  a 
NextEngine  HD  3-D  scanner.  Scanning  creates  a  replication  of  each  fragment,  which  preserves 
its  appearance  and  3-D  shape.  Scanning  must  be  completed  prior  to  any  required  cutting  for 
elemental  analysis.  The  fragment  is  placed  on  the  scanner  in  one  orientation,  is  scanned,  and 
then  is  rotated  into  an  orthogonal  position  in  relation  to  the  first  scan.  Then,  the  fragment  is 
scanned  again.  These  two  scans  are  then  merged  to  create  a  complete  360°  3-D  representation  of 
the  fragment  using  the  Geomagic  computer  program.  Any  holes  within  the  scans  due  to  defects 
in  the  exact  placement  of  the  fragment  or  in  the  processing  of  the  fragment  by  the  scanner  can  be 
filled  in  using  the  Geomagic  program.  Once  the  model  has  been  created,  “snapshots”  of  the 
model  are  taken  in  Geomagic  from  various  orientations  of  the  model  (top,  bottom,  left,  right, 
front,  back,  and  isometric).  The  final  snapshot  of  each  fragment,  the  file  containing  the  scan 
itself,  and  the  density  of  each  fragment  are  entered  into  the  fragment  database  (2). 

The  final  part  of  the  analysis  involves  the  elemental  analysis  of  each  fragment  by  scanning 
electron  microscopy-energy  dispersive  x-ray  (SEM-EDS),  which  may  be  accompanied  by 
inductively  coupled  plasma-atomic  emission  spectroscopy  (ICP-AES)  when  necessary.  A 
Hitachi  S-4700  field  emission  SEM  with  energy  dispersive  x-ray  spectroscopy  (EDS) 
capabilities  is  used  for  fragment  analysis.  A  low  magnification  (35x)  micrograph  is  taken  of  a 
representative  area  on  all  samples.  Qualitative  elemental  analysis  is  performed  with  EDS.  With 
this  technique  the  sample’s  surface  is  impacted  by  an  electron  beam.  The  electron  beam  excites 
the  surface  atoms  in  the  sample  producing  x-rays  that  are  characteristic  of  the  elements  found  in 
the  sample.  A  detector  is  used  to  convert  the  x-ray  energies  into  voltage  signals,  which  are  sent 
to  a  pulse-processor  for  detection.  The  results  are  displayed  as  a  spectrum  of  elements.  When 
indicated,  ICP-AES  may  be  used  to  determine  trace  concentrations  of  metal  and  is  used  to 
provide  quantitative  elemental  analysis  of  the  fragments. 

After  the  completion  of  all  analyses,  the  cases  are  transported  to  Range  10  on  Spesutie  Island  to 
be  stored  in  a  secure  evidence  cage.  In  the  past  two  months,  47  fragments  from  8  cases  were 
delivered  to  WSB  from  OAFME  and  analyzed  by  us  using  the  aforementioned  procedure. 


3.  Results  and  Discussion 


An  example  of  the  results  from  one  of  the  47  fragments  processed  from  the  past  two  deliveries  is 
illustrated  in  figure  1 .  The  physical  characteristics  of  the  fragment,  including  mass,  dimensions, 
recovery  location,  and  description,  were  taken  and  recorded  in  the  database.  Use  of  the  light 
microscope  revealed  multiple  tan  fibers  and  a  small  blue  fiber  adhering  to  the  surface  of  the 
fragment.  The  fragment  was  determined  to  be  1.33  g  with  dimensions  of  1 3 .93  x  14.42x3 .0 1  mm. 
Figure  2  provides  an  example  of  the  fragment  photograph  taken  with  a  forensic  evidence  L- 
shaped  ruler.  The  3-D  scan,  illustrated  in  figure  3,  preserved  the  image  and  shape  of  the 


101 


fragment  and  also  provided  the  fragment  density.  Scanning  indicated  the  density  was  2.20  g/mL. 
The  elemental  composition  of  this  fragment  was  determined  to  be  predominately  copper  and 
iron,  which  can  be  seen  in  the  SEM-EDS  spectrum  in  figure  4.  A  magnified  (35x)  surface  view 
of  an  area  on  the  fragment  is  shown  in  the  SEM  micrograph  (figure  5). 


Mass:  1.33  g 

Dimension:  13.93  x  14.42  x  3.01  mm 
Density:  2.20  g/mL 
Shape:  Irregular 

Recovery  Location:  Lateral  thigh 


Description:  Silver  colored  with  reddish- 
brown  encrusted  matter  and  ridges.  One 
end  is  more  jagged  and  irregular  whereas 
other  Is  more  square  with  rounded  edges. 
Many  tan  and  reddish  libers  on  convex  side 
along  with  a  small  blue  fiber  near  a  email 
hole  (located  near  area  where  fragment 
folded  onto  Itself)-  Near  Jagged  end  on 
reverse  side  are  several  tan  fibers  and 
fabric  with  spherical  adhesions  along 
concave  side. 

Predominant  Materials:  Copper  and  Iron 


Figure  1.  An  example  of  the  basic  information  entered  into  the  database.  The  density,  dimensions,  mass, 
shape,  recovery  location,  description,  and  predominant  materials  are  shown  for  this  fragment. 


Figure  2.  A  photographs  taken  of  the  fragment  once  sterilization 
has  been  completed. 
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Figure  3.  A  3-D  model  of  the  fragment  obtained  from 
the  scans  and  merged  using  Geomagic. 


Figure  4.  The  EDS  spectrum  from  the  qualitative  analysis  of  the  fragment  by  SEM-EDS. 
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Figure  5.  A  micrograph  of  the  fragment  obtained  by  SEM. 


4.  Summary  and  Conclusions 


The  elemental  composition  of  the  fragment  described  in  section  3  was  determined  to  be 
predominately  copper  and  iron.  These  results  were  combined  with  event  information  to 
determine  fragment  origin.  Fibers  were  found  on  the  fragment. 

The  results  of  the  analysis  of  these  47  fragments  were  provided  to  JTAPIC  partners  and  were 
briefed  to  JTAPIC’s  Dismounted  Soldier  Incident  Analysis  Team  (DIAT)  to  aid  in  their  event 
analysis.  The  3-D  models  and  specific  details  of  each  fragment  contained  in  the  JTAPIC 
database  will  also  be  used  in  future  modeling  and  simulation  efforts  and  in  event  recreations. 
WSB  analysts  will  also  use  the  results  to  increase  their  understanding  of  enemy  TTPs. 
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Abstract 


MBDyn  is  an  open  source  software  tool  capable  of  multi-body  dynamics  analysis.  The  Vehicles 
Technology  Directorate  (VTD)  of  the  U.S.  Army  Research  Laboratory  (ARL)  views  MBDyn  as 
a  potential  computational  research  and  design  tool  for  numerous  aerial  vehicle  platforms.  The 
most  common  application  of  MBDyn  is  for  rotor-blade  analysis,  a  subject  of  great  interest  to 
VTD.  However,  MBDyn  provides  a  very  generic  environment  for  creating  models  and 
conducting  dynamic  analysis,  making  it  suitable  for  numerous  applications.  One  such 
application,  also  of  interest  to  VTD,  is  research  concerning  flapping  wing  micro-aerial  vehicles 
(FWMAVs).  These  systems  are  driven  by  the  most  advanced  technologies  and  have  great 
potential  to  aid  the  operational  warfighter.  The  primary  objective  of  this  research  is  to  build  a 
preliminary  model  of  a  flapping  wing  system  in  MBDyn  and  evaluate  it  as  a  potential  tool  for 
future  use  at  VTD.  The  preliminary  models  of  flapping  wing  systems  presented  in  this  report 
illustrate  that  MBDyn  possesses  analysis  capabilities  well-suited  to  meet  many  VTD  needs. 
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1.  Introduction/Background 


Flapping  wing  micro-aerial  vehicles  (FWMAVs)  are  biologically  inspired  aircraft  micro-systems 
with  enormous  potential  for  Army  applications.  The  Vehicles  Technology  Directorate  (VTD)  of 
the  U.S.  Army  Research  Laboratory  (ARL)  in  recent  years  has  made  the  advancement  of  micro¬ 
systems  technology  one  of  their  objectives.  These  systems  would  be  extremely  useful,  discrete, 
and  efficient  surveillance  platforms,  which  would  enhance  the  warfighter’s  capabilities  and 
safety  while  in  operation.  Furthermore,  warfighter-focused  micro-systems,  especially 
FWMAVs,  are  driven  by  the  most  advanced,  cutting-edge  technologies,  which  bring  about  many 
design  challenges,  some  of  which  are  discussed  in  this  report.  The  most  obvious  and  critical 
design  challenge  for  FWMAVs  is  also  one  of  its  greatest  advantages:  there  are  seemingly  endless 
possibilities  for  FWMAV  configurations. 

The  vast  possibilities  are  easily  illustrated  by  noting  the  immense  variety  of  flapping  wing 
configurations  among  insects,  birds,  and  bats.  Figures  1  through  4  help  to  illustrate  the  variety 
and  complexity  of  wing  topologies.  Flapping  wing  flight  is  achieved  by  insects  with  masses  on 
the  order  of  milligrams  and  large  birds  with  masses  on  the  order  of  kilograms.  Furthermore, 
research  and  observation  has  shown  that  all  flapping  wing  species  achieve  flight  through 
different  patterns  of  wing  kinematics  and  different  wing  topologies.  There  are  a  vast  number  of 
variables,  all  of  which  play  some  role  in  determining  their  wing  shapes  and  wing  kinematics. 

For  example,  house  flies  depend  largely  on  forward  flight  as  opposed  to  hovering;  therefore,  the 
wings  of  a  house  fly  have  a  large  degree  of  sweep,  providing  increased  forward  thrust. 
Dragonflies  depend  on  forward  flight  as  well,  flying  several  miles  a  day,  but  also  need  the  ability 
to  hover;  therefore,  they  have  adapted  to  using  four  wings,  or  coupled  wings,  which  provide 
thrust  vectoring  capabilities  for  forward  flight,  hovering,  and  even  gliding.  Hummingbirds  are 
best  known  for  their  figure-of-eight  wing  kinematics  with  very  high  wing  beat  frequencies.  This 
design  provides  them  with  enhanced  maneuverability  and  agility.  Most  other  birds  have  wing 
kinematics  with  a  much  more  limited  range  of  motion,  and  typically  have  the  lowest  wing  beat 
frequencies  amongst  flapping  wing  species.  Another  reason  for  the  immense  variety  is  that  the 
aerodynamic  effects  change  drastically  over  a  broad  range  of  low  Reynolds’  numbers. 

Thousands  of  years  of  natural  selection  have  fine-tuned  these  variables,  endowing  each  particular 
species  with  the  flight  capabilities  necessary  for  survival.  For  successful  development  of 
FWMAV  technology,  the  ultimate  goal  is  to  understand  the  sensitivities  that  various  flight 
objectives  have  toward  all  of  these  variables. 
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Figure  1.  Dragonfly  (7). 


Figure  2.  Wing  topology  (2). 


Figure  3.  Wing  nomenclature  (3). 


Figure  4.  Dragonfly  wing  venation  ( 4 ). 
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At  the  prime  of  the  aviation  industry,  engineers  and  inventors  approaching  the  problem  of  flight 
were  faced  with  similar  challenges.  In  order  to  design  a  successful  aircraft  system,  they  needed 
to  understand  the  numerous  physical  relationships  in  play.  After  a  century’s  worth  of  knowledge 
and  experience,  these  relationships  for  conventional  aircraft  systems  are  somewhat  trivial  and  the 
design  process  has  been  drastically  simplified  to  a  fundamental  science.  With  interest  growing 
in  unconventional  aircraft  systems,  such  as  FWMAVs,  all  of  the  previous  knowledge  and 
experience  cannot  be  directly  applied.  Though  it  can  be  very  helpful,  the  result  is,  once  again, 
we  face  a  very  complicated  and  uncharted  design  process. 

Experimental  studies  will  be  very  important  to  designing  successful  flapping  wing  (FW) 
platforms,  but  computational  studies  will  be  crucial  for  optimization.  Computational  analysis  of 
FW  flight  depends  heavily  on  dynamics  and  aeroelastic  effects.  For  these  systems,  optimization 
does  not  come  in  the  form  of  mass  minimization  and  stress  constraints,  like  it  does  for 
conventional  aircraft  systems,,  rather,  optimization  of  FWMAVs  depends  on  aeroelastic  tailoring 
of  the  wing  topology  to  optimize  certain  flight  parameters,  such  as  cruise  efficiency  or  hover 
capability.  The  optimization  can  be  constrained  by  controllability  metrics,  manufacturing 
constraints,  etc.  Therefore,  it  is  imperative  that  an  efficient  and  capable  computational  design 
tool  be  identified.  The  primary  objective  of  this  research  is  to  evaluate  a  computational  design 
tool  called  MBDyn.  This  multi-body  dynamics  tool  has  the  potential  to  be  used  for  detailed 
design  of  FWMAVs  and  other  aerial  vehicle  platforms  across  the  board.  Currently,  there  are 
several  documented  applications  of  rotor  blade  analysis  in  MBDyn,  all  of  which  speak  well  for 
MBDyn’s  capabilities.  However,  there  is  no  documentation  of  MBDyn  being  used  to  model  FW 
systems.  Therefore,  this  research  evaluates  the  efficiencies  and  inefficiencies  of  MBDyn  through 
a  few  simple  FW  models. 


2.  Model  Definition 


VTD  is  developing  experimental  models  of  several  micro-system  platforms  it  is  researching.  Dr. 
Asha  Hall  of  the  Multifunctional  Structures  team  in  the  Mechanics  Division  has  built  an 
experimental  FW  model  that  makes  use  of  piezoelectric  (lead  zirconate  titanate  [PZT])  actuation. 
Hall’s  model  is  used  as  the  baseline  for  developing  a  FW  model  in  MBDyn.  The  following 
sections  describe  both  the  experimental  model  and  the  computational  model,  while  also  detailing 
MBDyn. 

2.1  Experimental  Flapping  Wing  Model 

PZTs  are  often  used  for  vibration  damping  on  a  variety  of  rotor-craft  platforms  and  in  many 
other  applications  as  well.  In  this  case,  PZT  is  being  used  for  the  exact  opposite  reason,  to  excite 
motion  in  the  system.  The  PZT  in  this  case,  is  a  thin  cantilevered  beam  and  can  be  seen  clamped 
in  the  orange  vice  in  figure  5.  One  end  of  the  beam  is  clamped  and  connected  to  a  power  source. 
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The  beam  has  three  layers  of  material,  and  as  voltage  is  applied,  the  top  and  bottom  layer  act  as 
opposite  poles,  while  the  middle  layer  is  a  conductive  passageway  for  electron  transfer.  The 
polarity  between  the  two  beams  creates  a  strain,  causing  the  PZT  to  deflect.  If  the  poles  are 
reversed,  the  direction  of  the  strain  reverses,  thus  deflecting  the  beam  in  the  opposite  direction. 
Therefore,  the  end  of  the  beam  can  be  forced  to  oscillate  up  and  down  at  a  desired  frequency. 
The  particular  model  that  Hall  has  built  has  two  wings  rooted  at  the  oscillatory  end  of  the  PZT. 
Thus,  the  up  and  down  oscillation  of  the  PZT  excites  a  flapping  motion  in  the  wings.  This  is 
becoming  a  very  popular  method  of  actuation  for  FW  systems  due  to  its  high  frequency  response 
and  the  ability  to  map  periodic  oscillations  into  desired  wing  kinematics  (5,  6). 


Figure  5.  MRI  experimental  FW  model  (5). 

Hall  created  the  wing  structure  using  a  three-dimensional  (3-D)  printing  method.  The  structure 
was  modeled  in  computer-aided  design  (CAD)  software,  in  this  case,  SolidWorks.  The  CAD 
model  was  then  provided  to  the  3-D  printer,  in  which  a  FullCure  cartridge  was  installed. 

FullCure  was  the  material  used  to  create  the  structure  and  it  is  commonly  used  for  3-D  printing. 
A  couple  hours  later,  Hall  was  able  to  remove  her  FullCure  structure  from  the  printer  and  begin 
conducting  her  experiments.  The  wing  surface  is  a  low  density  polyethylene  (LDPE)  film.  In 
figure  5,  note  the  gray  patches  near  the  wing  tips.  These  are  small  pieces  of  reflective  tape, 
which  give  Hall  the  ability  to  take  displacement  measurements  while  the  wings  are  flapping  (5). 

Currently,  the  wings  are  connected  via  the  two  carbon  fiber  rods  seen  in  figure  5.  These  rods  are 
connected  in  such  a  way  that  they  oscillate  in  unison  with  the  PZT.  The  wing  connection 
behaves  like  a  clamped  beam  and  does  not  allow  rotation.  The  wings  are  flexible  enough  that  if 
excited  at  the  right  frequency  they  will  take  on  a  flapping  motion.  Future  models  that  MRI 
builds  will  look  towards  creating  a  mechanical  connection  that  allows  rotation.  This  will  allow 
the  flapping  motion  to  achieve  larger  amplitudes  than  the  current  model.  In  addition,  the  current 
model  only  gives  Hall  active  control  of  wing  plunge.  Many  research  studies  have  shown  that  the 
phase  difference  between  wing  plunge  and  wing  pitch  is  key  to  optimizing  thrust  production. 
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Therefore,  MRI  is  in  the  process  of  designing  a  modified  PZT  bimorph  that  will  give  them  active 
control  of  both  wing  plunge  and  wing  pitch.  This  will  allow  MRI  to  experimentally  explore  a 
broader  range  of  wing  kinematics  (5). 

2.2  MBDYN  Input 

MBDyn  runs  off  of  a  text  input  file  that  requires  specific  formatting.  The  input  file  is  divided 
into  five  different  data  sets.  This  section  briefly  defines  those  data  sets  and  the  basic  layout  for 
MBDyn  input.  A  partial  example  of  an  MBDyn  input  file  is  provided  in  appendix  A.  Section 

2.3  provides  the  detail  of  how  the  FW  model  specific  to  this  research  was  defined  in  MBDyn. 

The  five  data  sets  are  problem  type,  problem,  control  data,  nodes,  and  elements.  (7) 

The  first  data  set  simply  tells  MBDyn  what  type  of  problem  it  will  be  solving.  The  analysis  used 
for  the  FW  model  was  an  initial  value  problem;  therefore,  the  second  data  set,  specific  to  the 
problem  type,  is  called  initial  value.  This  data  set  requires  the  definition  of  several  parameters 
and  allows  for  the  optional  definition  of  several  others.  The  required  parameters  are  initial  time, 
final  time,  time  step,  max  iterations,  and  tolerance.  The  initial  and  final  time  simply  define  the 
time  domain  within  which  the  analysis  is  conducted.  The  time  step  is  the  increment  size  by 
which  the  time  domain  is  discretized.  Max  iterations  limits  the  convergence  of  the  solution  at 
each  time  step  to  a  maximum  number  of  iterations.  Tolerance  is  the  absolute  error  that  must  be 
achieved  at  each  time  step.  (7) 

The  third  data  set  is  called  control  data.  This  set  informs  MBDyn  of  what  is  included  in  the 
model  definition.  That  is,  what  type  and  how  many  nodes  and  elements  are  used  in  the  model. 
This  prepares  MBDyn  for  what  to  look  for  and  read  and  may  not  be  required  in  future  versions  of 
MBDyn,  but  it  certainly  helps  for  readability  and  debugging  purposes.  MBDyn’s  definition  of 
node  and  element  covers  many  types  of  nodes  and  elements.  For  instance,  a  node  can  be 
structural,  abstract,  hydraulic,  etc.,  and  an  element  is  considered  to  be  a  beam,  a  joint,  a  rigid 
body,  a  force,  etc.  (7) 

The  fourth  data  set  is  called  nodes.  All  MBDyn  nodes  used  in  the  model  are  defined  in  this  data 
set.  For  the  FW  model,  all  nodes  are  defined  as  structural  nodes;  however,  there  are  a  few 
different  types  of  structural  nodes,  most  notably,  static  and  dynamic.  Basically,  a  dynamic  node 
can  be  issued  mass  whereas  a  static  node  cannot.  In  the  FW  model,  there  is  one  static  node  at  the 
wing  root,  necessary  to  ground  the  entire  structure,  while  all  other  nodes  are  dynamic.  The 
definition  of  a  node  requires  the  following  information:  node  type,  node  ID,  absolute  position, 
absolute  orientation,  absolute  velocity,  and  absolute  angular  velocity.  Here  the  word  “absolute” 
refers  to  the  global  reference  frame  or  the  absolute  reference  frame.  (7) 

The  fifth  data  set  is  called  elements.  All  MBDyn  elements  used  in  the  model  are  defined  in  this 
data  set.  For  the  FW  model,  there  were  four  different  types  of  elements  used:  joints,  rigid 
bodies,  beams,  and  forces.  MBDyn  has  a  joint  library,  which  allows  one  to  define  everything 
from  a  clamp  joint  to  a  pin  joint  to  a  revolute  hinge  to  a  gimbal  hinge.  This  library  is  diverse 
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enough  that  almost  any  imaginable  mechanical  connection  can  be  modeled.  This  versatility  will 
be  important  for  future  development  of  FW  models  with  complicated  actuation  and  wing 
kinematics.  This  particular  FW  model  makes  use  of  a  clamp  joint,  connecting  the  wing  root  to 
the  leading  edge  spar,  and  a  total  joint.  The  total  joint  provides  the  ability  to  constrain  or  drive 
any  degrees  of  freedom  of  the  relative  node.  A  drive  can  be  defined  as  any  time-dependent 
function  and  prescribed  to  a  specific  degree  of  freedom  of  any  node.  The  drive  used  to  prescribe 
motion  in  the  FW  model  is  described  in  section  2.3.  Rigid  body  definition  gives  the  structure 
mass  and  inertial  properties.  Any  number  of  rigid  bodies  can  be  defined,  and  their  definition 
consists  of  a  rigid  body  ID,  the  node  ID  relative  to  which  the  rigid  body  is  being  defined,  the 
total  mass  of  the  rigid  body,  a  position  vector  defining  the  center  of  mass  of  the  rigid  body,  and  a 
mass  moment  of  inertia  matrix  for  the  rigid  body.  This  gives  users  versatility  in  the  type  of  rigid 
body  they  wish  to  model.  Beams  are  an  element  type  that  defines  nodal  connectivity.  Beam2  is 
a  beam  connecting  two  nodes,  while  beam3  is  a  beam  connecting  three  nodes.  The  Beam2 
definition  requires  a  beam  ID;  the  two  node  IDs  which  the  beam  is  connecting,  a  relative 
orientation  matrix,  which  allows  the  beam  to  be  defined  in  a  local  reference  frame;  and  a 
constitutive  law  by  which  the  beam  behaves.  (7) 

2.3  MBDYN  Flapping  Wing  Model 

The  first  couple  of  weeks  of  this  research  time  were  spent  climbing  the  initial  learning  curve  of 
MBDyn.  Very  simple  cantilever  models  were  built  and  analyzed  in  MBDyn,  and  the  research 
slowly  progressed  to  the  point  of  modeling  a  FW  system.  In  the  end,  a  rather  refined  process 
was  developed  to  build  a  FW  model  for  analysis  in  MBDyn. 

The  process  begins  by  creating  a  two-dimensional  unstructured  surface  mesh  of  triangular 
elements  using  a  program  called  Solid  Mesh.  For  the  baseline  model  (MRI’s  experimental 
model),  the  unstructured  mesh  was  created  by  mapping  it  over  the  SolidWorks  file  that  Hall  used 
for  printing  the  FullCure  structure.  This  ensured  that  the  dimensions  of  the  wing  boundaries  for 
the  computational  model  matched  those  of  the  experimental  model.  Figures  6  and  7  show  the 
SolidWorks  model  of  the  FullCure  structure  and  the  resulting  unstructured  mesh. 


Figure  6.  Image  of  the  SolidWorks  wing  model  (5). 
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Figure  7.  Unstructured  mesh  for  the  MBDyn  model  construction. 

From  here,  the  user  employs  a  pool  of  MATLAB  scripts  that  were  developed  throughout  this 
research.  The  first  MATLAB  script  reads  the  Solid  Mesh  data  file  and  stores  the  element  data 
into  the  workspace.  In  addition  to  the  element  data,  some  user  input  is  required.  The  user  is 
prompted  to  provide  the  following  information: 

1 .  The  Node  ID  number  at  which  the  hinge  is  to  be  located. 

2.  The  number  of  material  groups  the  user  wishes  to  define. 

3.  Material  properties  for  the  default  material  group  (prescribed  to  all  structural  members): 

a.  Material  group  ID 

b.  Elastic  modulus 

c.  Shear  modulus 

d.  Mass  density 

e.  Thicknesses  of  member  cross  sections 

4.  Properties  1-4  for  each  subsequent  material  group 

5.  The  number  of  stiffening  members  the  user  wishes  to  define  (spars/ribs). 

6.  For  each  stiffening  member,  the  user  must  enter  the  following: 

a.  Material  group  ID  used  for  the  stiffening  member 

b.  Vector  of  nodal  connectivity  for  the  stiffening  member 

c.  Thicknesses  of  the  member  cross  section  at  first  node 

d.  Thicknesses  of  the  member  cross  section  at  last  node 

The  material  property  for  this  particular  FW  model  is  the  polyethylene  material  group.  These 
material  properties  are  prescribed  to  all  structural  members.  The  material  properties  list  in  item  4 
are  used  to  define  the  FullCure.  There  are  two  stiffening  members  defined  in  this  FW  model: 
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one  is  the  leading  edge  FullCure  spar  and  the  other  is  the  root  edge  FullCure  rib.  A  linear 
interpolation  is  applied  to  the  thicknesses  defined  in  item  6,  allowing  the  leading  edge  spar  and 
root  edge  to  be  properly  tapered. 

A  second  MATLAB  script  takes  the  user  input  and  element  data  and  writes  all  of  the  information 
to  an  MBDyn  input  file  with  proper  formatting.  The  MATLAB  script  then  calls  and  executes 
MBDyn.  If  MBDyn  has  an  error,  it  prints  the  error  to  the  MATLAB  command  window,  and  if  it 
is  successful,  another  MATLAB  script  reads  the  output  files  and  stores  the  information  into  the 
MATLAB  workspace.  The  MATLAB  script  has  some  optional  post-processing  features  for 
graphs  and  visualization.  In  addition  to  the  MATLAB  visualization  shown  in  figure  8,  there  is 
an  application  called  EasyAnim,  which  is  compatible  with  MBDyn.  A  snapshot  of  the 
EasyAnim  window  is  provided  in  figure  9.  EasyAnim  allows  for  the  user  to  create  video  files  to 
visualize  the  structural  dynamics  for  the  entire  time  domain  of  the  analysis.  Before  presenting 
the  results  of  this  model,  it  is  important  to  discuss  some  modeling  considerations. 


Figure  8.  MATLAB  visualization  example. 
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Figure  9.  EasyAnim  visualization  snapshot. 


2.4  Important  Modeling  Considerations 

This  section  discusses  important  issues  regarding  beam  analysis,  equivalent  continuum,  rigid 
body  definition,  and  wing  excitation.  All  of  these  topics  have  considerable  influence  on  the 
dynamic  response  of  the  system. 


As  seen  in  the  figure  9,  the  wing  surface  is  being  modeled  as  a  network  of  one-dimensional  (1-D) 
beams.  Before  discussing  the  issue  of  equivalent  continuum,  it  is  important  to  understand  how 
MBDyn  is  conducting  the  deformation  analysis  of  1-D  beams.  For  the  FW  application,  like 
many  others,  being  able  to  accurately  model  large  deformations  and  rotations  is  critical.  It  is 
here  that  many  beam  theories  can  lose  accuracy.  The  beam  theory  being  used  by  MBDyn  is  the 
geometrically  exact  beam  theory  (GEBT).  GEBT  is  essentially  the  combination  of  a  1-D 
nonlinear  beam  analysis  and  a  2-D  cross-sectional  stiffness  analysis.  It  relies  on  the  constitutive 
law  (equation  1);  however,  the  constant  matrix  does  not  have  to  be  diagonal  and  can  include  the 
proper  coupling  terms.  In  addition,  variations  of  this  law  can  be  used  that  include  warping  and 
viscous  effects  ( 8 ,  9). 
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This  law  allows  for  direct  computation  of  the  strains.  The  strains  are  then  simply  integrated  in 
order  to  calculate  the  displacements  and  rotations.  For  that  reason,  this  type  of  beam  analysis 
can  capture  large  deformations  while  only  relying  on  the  assumption  of  small  strains.  The  FW 
system  is  perfect  for  this,  because  flexible  wings  will  respond  with  large  rotations  in  the  local 
reference  frame  but  very  small  strains.  Furthermore,  the  constitutive  law  contains  the  entire 
geometric  definition  of  the  beam,  making  no  approximations,  hence  the  name  “geometrically 
exact  (8,  9).” 

For  the  purposes  of  this  research,  it  would  be  ideal  to  represent  the  wing  membrane  with  either  2- 
D  shell  elements  or  membrane  elements.  Unfortunately,  as  previously  stated,  MBDyn  is  limited 
to  1-D  beams.  So,  it  is  important  to  consider  how  a  beam/truss  lattice  continuum  can  be  equated 
to  a  thin  film/membrane  continuum.  Odegard  and  Gates  (10)  present  an  equivalent  continuum 
method  for  representing  a  complicated  volume  with  a  finite  element  truss  lattice.  In  their  paper, 
Odegard  and  Gates  define  certain  criteria,  which  if  met,  means  the  representative  truss  lattice  can 
be  assumed  equivalent  to  the  desired  continuum.  When  the  two  models  are  prescribed,  a  discrete 
load  under  static  analysis  must  satisfy  the  following: 

•  The  boundary  displacement  of  the  representative  element  must  match  that  of  the  other 
model. 


•  The  thermo-elastic  strain  energies  of  the  two  models  must  be  equivalent. 

Odegard  and  Gates  provide  the  following  expression  for  calculating  the  strain  energy  of  a 
discrete  beam/truss  lattice. 


A1 


(2) 


Here,  A  is  the  cross-sectional  area  of  the  member,  Y  is  Young’s  modulus  of  the  member  material, 
r  is  the  deformed  length  of  the  member,  and  R  is  the  un-deformed  length  of  the  member.  The 
summation  about  j  is  from  1  to  the  number  of  types  of  members,  if  there  are  multiple  types;  and 
the  summation  about  i  is  from  1  to  the  number  of  members  (10). 


The  primary  issue  with  applying  this  method  to  the  model  defined  in  this  research  is  that  the 
method  depends  upon  periodicity  of  the  representative  lattice.  Unfortunately,  the  unstructured 
mesh  used  to  create  this  model  does  not  have  periodicity.  Nevertheless,  a  lot  can  be  said  for  the 
importance  of  equating  the  strain  energies  of  two  different  models  in  order  to  consider  them 
equivalent  continua.  The  results  presented  in  the  research  do  not  attempt  to  apply  an  equivalent 
continuum  method  to  the  model.  Rather,  this  preliminary  computational  model  uses  the  material 
properties  of  polyethylene  while  keeping  the  representative  beam  members  thin  in  cross  section. 
Equating  the  strain  energy  of  the  representative  model  to  the  polyethylene  membrane  should  be 
applied  in  the  future,  and  is  further  discussed  in  section  4  of  this  report  (10-12). 
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Rigid  body  definition  could  be  done  in  a  number  of  different  ways  and  still  properly  represent 
the  network  of  1-D  beams;  however,  the  simplest  way  is  to  define  a  rigid  body  for  each  beam. 
Thus,  when  each  beam  is  identified  with  two  nodes,  the  rigid  body  for  that  beam  is  defined 
relative  to  the  first  of  the  two  nodes.  The  mass  is  defined  by  multiplying  the  mass  density  of  the 
material  for  each  beam  by  its  respective  length  and  cross-sectional  area.  The  relative  center  of 
mass  is  defined  as  the  distance  vector  from  the  relative  node  to  the  midpoint  of  the  beam.  Lastly, 
the  mass  moment  of  inertia  matrix  is  defined  as  that  of  a  simple  rod  or  beam  in  the  absolute 
reference  frame  (9). 

The  last  key  component  of  the  model  definition  is  prescribing  the  dynamic  motion  of  the  system. 
All  nodes  are  initialized  at  rest,  and  the  wing  surface  is  defined  in  the  x-y  plane  with  zero  pitch. 
As  previously  described,  the  current  wing  root  connections  force  the  wings  to  have  the  behavior 
of  clamped  boundary  conditions.  This  dynamic  motion  is  easily  modeled  in  MBDyn  through 
what  is  called  a  drive.  A  drive  is  a  time-dependent  function,  usually  periodic,  which  defines  the 
motion  of  a  given  node.  Consequently,  in  this  FW  model,  the  drive  was  defined  as  a  sine 
function,  which  acted  for  the  entire  time  domain.  As  previously  mentioned,  the  current 
experimental  model  has  wings  that  behave  as  clamped  beams  constraining  rotation  at  the  roots. 
Therefore,  two  of  the  models  presented  in  section  3  excited  a  vertical  oscillation  of  the  wing  root. 
In  addition,  section  3  also  presents  a  model  that  was  excited  by  prescribing  a  rotation  at  the  root 
rather  than  a  vertical  oscillation. 


3.  Results 


This  section  details  three  very  simple  FW  models.  These  models  illustrate  that  MBDyn  is 
capable  of  successfully  modeling  a  simple  FW  system.  Section  4  of  this  report  sheds  light  on 
MBDyn’ s  potential  role  for  conducting  computational  studies  to  benefit  and  advance  the  efforts 
of  detailed  design  of  FWMAV  platforms.  No  comparison  of  the  computational  results  to  the 
experimental  results  is  made  here.  Some  refinement  needs  to  be  done  before  the  computational 
model  can  be  directly  compared  to  the  experimental  model.  Right  now,  the  wing  root  excitation 
of  the  computational  model  is  of  arbitrary  amplitude.  This  is  because  displacement 
measurements  of  the  PZT  have  not  been  directly  taken.  Thus,  the  dynamic  motion  prescribed  in 
the  computational  model  cannot  yet  be  equated  to  that  of  the  experimental  model.  Nonetheless, 
the  results  presented  illustrate  MBDyn’s  capability  to  model  FW  systems. 

3.1  Model  1 

The  first  model  is  built  upon  the  coarse,  unstructured  mesh  shown  in  figure  10.  This  mesh  was 
mapped  directly  from  the  SolidWorks  model  used  to  print  the  FullCure  structure  of  the 
experimental  model. 
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Figure  10.  Coarse,  unstructured  mesh. 


The  following  are  the  specifications  for  the  model: 


Membrane  material: 

Elastic  modulus: 

150  MPa 

Shear  modulus: 

150  MPa 

Mass  density: 

1030  kg/mA3 

Cross  section: 

1  mm  x  1  mm 

FullCure  material: 

Elastic  modulus: 

700  MPa 

Shear  modulus: 

700  MPa 

Mass  density: 

1092  kg/mA3 

Leading  edge  spar: 

Material: 

FullCure 

Cross  section  at  node  15: 

1 .29  mm  x  1 .26  mm 

Cross  section  at  node  1 : 

1  mm  x  0.52  mm 

Root  Edge: 

Material: 

FullCure 

Cross  section  at  node  7: 

1  mm  x  0.66  mm 

Cross  section  at  node  15: 

1 .29  mm  x  1 .26  mm 

Before  exciting  the  structure  and  conducting  a  dynamic  analysis,  it  is  important  to  conduct  an 
eigenvalue  analysis  of  the  structure.  Every  structure  has  various  natural  frequencies.  When  a 
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structure  is  excited  at  one  of  its  natural  frequencies  it  will  respond  according  to  the 
corresponding  mode  shape.  In  this  particular  case,  a  flapping  motion  is  desired,  and  twist  about 
the  span-wise  axis  is  not  desired.  After  conducting  an  eigenvalue  analysis  of  this  structure, 
another  capability  of  MBDyn,  the  mode  shapes  shown  in  figure  1 1  resulted. 


19.8  Hz  71.3  Hz 

Figure  11.  Model  1  mode  shapes. 
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111.7  Hz 
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The  lowest  frequency  results  in  the  lowest  degree  of  curvature  with  little  twist.  The  higher 
frequencies  exhibit  motion  that  has  either  too  much  twist  or  too  much  curvature,  neither  of  which 
is  desired.  Therefore,  when  conducting  this  dynamic  analysis,  the  excitation  should  be  at  the 
first  frequency. 


Figures  12  through  15  illustrate  the  structural  response  of  the  wing  due  to  an  excitation  at  a 
frequency  of  19.8  Hz.  The  excitation  was  prescribed  at  the  wing  root  in  the  vertical  z-direction 
with  an  amplitude  of  1  mm.  Notice  that  the  wing  tip  displacement  peaks  at  10  mm.  Because  the 
excitation  was  at  a  natural  frequency,  the  response  is  amplified  as  the  excitation  continues, 
inducing  a  flapping  motion  in  the  wing.  It  should  also  be  noted  that  this  model  has  zero 
damping,  which  will  be  necessary  for  future  models.  In  addition,  a  beat  frequency  is  apparent, 
which  is  the  result  of  the  interaction  between  multiple  frequencies.  This  is  less  than  desirable, 
and  will  have  to  be  explored  in  future  models. 


Figure  12.  Model  1  excitation,  amplitude  =  1  mm,  frequency  =  19.8  Flz. 
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Figure  13.  Model  1  wing  tip  displacement. 


Figure  14.  Model  1  wing  tip  rotations. 


128 


1.5 

1 

\  0.5 

B 

£?  o 

u 

_o 

£  -0.5 

-1 

-1.5 

C 

Wing  Tip  Velocity 

A . L.  . 

A . A  A 

ft 

aJ  I 

A  A 

i 

III 

I  :  [  l  /  V  \  1  j: 

I  I  I/  "  1; 

II  II  V v 

1  f  '1 

J  i 

i 

1 

0.2  0.4  0.6  0.8  1 

Time,  s 

Figure  15.  Model  1  wing  tip  velocity. 

3.2  Model  2 

The  second  model  is  built  upon  the  fine,  unstructured  mesh  shown  in  figure  16.  This  mesh  was 
also  mapped  directly  from  the  SolidWorks  model,  but  given  a  smaller  element  spacing. 


Wing  Root:  Node  15 


Root  Edge 


Node  7 


Span  (m) 


Figure  16.  Fine,  unstructured  mesh. 

The  following  are  the  specifications  for  the  model: 
•  Membrane  material: 


Elastic  modulus: 
Shear  modulus: 
Mass  density: 
Cross  section: 


150  MPa 
150  MPa 
1030  kg/mA3 
1  mm  x  1  mm 


Leading  Edge  Spar 


Node  1 
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FullCure  material: 

Elastic  modulus: 

700  MPa 

Shear  modulus: 

700  MPa 

Mass  density: 

1092  kg/mA3 

Leading  edge  spar: 

Material: 

FullCure 

Node  connectivity: 

[15,31:53,1] 

Cross  section  at  node  15: 

1 .29  mm  x  1 .26  mm 

Cross  section  at  node  1 : 

1  mm  x  0.52  mm 

Root  Edge: 

Material: 

FullCure 

Node  connectivity: 

[7,22:30,15] 

Cross  section  at  node  7: 

1  mm  x  0.66  mm 

Cross  section  at  node  15: 

1 .29  mm  x  1 .26  mm 

The  eigenvalue  analysis  was  conducted  the  same  as  in  Model  1 ;  however,  the  results  this  time 
were  less  than  desirable.  As  can  be  seen  from  the  mode  shape  in  figure  17,  no  natural  frequency 
exists  that  produces  the  flapping  motion  desired.  This  model  was  tested  at  the  lowest  natural 
frequency  but  resulted  in  a  motion  involving  little  wing  plunge  and  a  lot  of  wing  twist/pitch. 


Figure  17.  Model  2  mode  shape. 

This  is  a  perfect  illustration  of  how  important  equivalent  continuum  is  for  continuing  research. 
Ultimately,  an  equivalent  continuum  method  would  allow  for  a  particular  system  to  be  accurately 
modeled  by  different  meshes.  The  finer  mesh  used  in  this  model  makes  the  wing  much  stiffer 
due  to  the  larger  number  of  beams.  A  stiffer  structure  results  in  higher  natural  frequencies,  and 
in  this  case,  the  mode  shapes  result  in  motions  that  are  not  desired.  Since  an  accurate  equivalent 
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continuum  method  has  yet  to  be  developed,  all  that  can  be  done  is  adjust  the  beam  properties  by 
way  of  trial  and  error  until  a  similar  frequency  and  mode  shape  is  obtained.  Upon  doing  so  the 
membrane  properties  were  adjusted  to  E  =  G  =  132  Mpa,  and  new  eigenvalues  and  mode  shapes 
were  found  (figure  18). 


18.5  Hz  56.0  Hz 


65.4  Hz  102.3  Hz 


172.8  Hz 


Figure  18.  Model  2  mode  shapes  with  adjusted  membrane  properties,  E  =  G  =  132  Mpa. 

Now  the  structure  has  a  natural  frequency  that  corresponds  to  a  mode  shape  very  similar  to  that 
of  Model  1 .  The  structure  was  excited  at  18.5  Hz  with  the  same  amplitude  used  in  Model  1 .  The 
results  are  shown  in  figures  19  through  22.  The  motion  exhibits  the  same  behavior  as  Model  1, 
but  notice  that  the  maximum  displacement  of  the  wing  tip  has  decreased  to  about  8.5  mm.  For 
future  models,  it  will  be  very  critical  to  represent  the  membrane  appropriately.  Changing  the 
mesh  forced  a  change  in  the  material  properties  to  obtain  similar  results.  These  sensitivities  will 
need  to  be  explored  in  addition  to  formulating  an  adequate  equivalent  continuum  method. 


Figure  19.  Model  2  excitation,  amplitude  =  1  mm,  frequency  =  18.5  Flz. 
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Figure  20.  Model  2  wing  tip  displacement. 
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Figure  21.  Model  2  wing  tip  rotations. 
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Figure  22.  Model  2  wing  tip  velocity. 

3.3  Model  3 

The  third  model  varies  from  the  first  two  in  that  the  motion  prescribed  at  the  wing  root  is  a 
rotational  motion  as  opposed  to  a  vertical  displacement.  This  excitation  forces  the  orientation  of 
the  wing  to  change  more  drastically  than  that  of  the  first  two  models.  This  model  was  built  upon 
the  coarse  mesh  used  in  Model  1  and  has  all  of  the  same  material  properties.  Because  it  was  the 
same  structure  as  Model  1,  it  has  the  same  natural  frequencies,  and  so,  was  excited  at  the  same 
frequency  of  19.8  Hz.  The  results  are  shown  in  figures  23  through  26.  Because  rotation  was 
allowed,  the  wing  tip  displacement  achieved  a  larger  amplitude,  exhibiting  a  more  enhanced 
flapping  motion. 


Figure  23.  Model  3  excitation,  amplitude  =  0.1  radians,  frequency  =  19.8  Flz. 
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Wing  Tip  Displacement 


Figure  24.  Model  3  wing  tip  displacement. 


Figure  25.  Model  3  wing  tip  rotations. 
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Wing  Tip  Velocity 
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Figure  26.  Model  3  wing  tip  velocity. 


4.  Discussion  and  Future  Work 


4.1  Short-term  Objectives 

The  primary  short-term  objective  is  to  obtain  the  displacement  data  for  the  PZT  in  the 
experimental  model.  This  exact  same  motion  could  then  be  applied  to  the  computational  model. 
Making  sure  that  all  other  material  parameters  were  consistent,  the  experimental  and 
computational  models  could  then  be  compared.  This  would  serve  as  a  validation  for  the  MBDyn 
model  as  well  as  lead  into  other  interesting  studies  involving  the  comparison  of  MBDyn  results 
and  experimental  results.  Models  to  be  explored  consist  of  varying  the  type  of  excitation 
function.  Due  to  the  beating  seen  in  these  models,  it  is  apparent  that  the  excitation  amplifies  the 
motion  to  a  certain  limit  and  then  it  interferes.  In  addition,  adding  pitch  to  the  mix  should  create 
for  some  interesting  results. 

As  pointed  out  several  times,  modeling  the  membrane  accurately  with  1-D  beams  is  going  to  be  a 
challenge.  Trade  studies  will  be  done  to  detennine  how  sensitive  the  response  is  to  mesh 
refinement,  changing  the  value  of  EA  that  is  prescribed  to  the  representative  beam  elements,  etc. 
This  could  consist  of  scaling  Young’s  modulus  as  well  as  the  cross-sectional  area  of  the 
representative  beams. 

Another  very  important  short-term  objective  is  to  begin  the  process  of  pairing  MBDyn  with  an 
external  computational  fluid  dynamics  (CFD)  tool.  An  MBDyn  model  can  be  set  up  to 
communicate  with  external  software  via  an  edge  communicator.  Essentially,  MBDyn  is  told  to 
print  specified  data  at  every  time  step  iteration  to  a  data  file.  The  external  CFD  tool  does  the 
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same,  and  the  two  software  share  this  file,  which  is  called  an  edge  communicator.  The 
communicator  tells  both  MBDyn  and  the  CFD  tool  when  to  read  data  and  when  to  stop  reading 
data.  The  long-temi  objectives  discussed  in  section  4.2  all  depend  on  this  capability.  Detailed 
design  of  FW  models  requires  accurate  aeroelastic  analysis  of  the  system  dynamics.  In  addition, 
the  aerodynamics  involved  with  FW  flight  include  3-D  unsteady  effects  dominated  by  vortex 
flows  and  wing-wake  interaction,  thus  requiring  an  analysis  capable  of  fluid-structure 
interaction.  This  means  conducting  both  structural  and  fluid  analysis  simultaneously  because  the 
output  of  one  analysis  is  the  input  to  the  other  and  vice-versa.  In  addition,  including  the 
aerodynamic  forces  will  act  as  a  system  damper.  This  will  play  a  key  role  in  having  the  system 
approach  a  limit  cycle. 

4.2  Long-term  Objectives 

The  ultimate  objective  is  to  a  have  methodical  approach  to  the  detailed  design  of  FWMAVs.  As 
previously  mentioned,  the  aeroelastic  effects  involved  with  flexible  wing  configurations  require 
an  analysis  capable  of  fluid-structure  interaction.  Once  MBDyn  is  coupled  with  an  external  CFD 
solver,  the  computational  design  tool  will  be  completely  in  place.  Optimization  is  the  next  big 
design  step,  but  in  order  to  do  this  effectively,  classifying  and  parameterizing  various  FW 
platforms  is  crucial.  Optimization  depends  on  objective  functions,  design  variables,  and 
constraints,  and  identifying  and  defining  such  variables  and  constraints  is  a  challenge. 

As  discussed  in  section  1,  different  FW  species  depend  on  different  flight  characteristics.  A  fly 
depends  on  quick  forward-flight,  a  hummingbird  depends  on  hovering  and  agility,  while  a 
bumblebee  depends  on  efficient  cruising  and  hovering.  Therefore,  different  FW  platforms  will 
depend  on  different  flight  functions  as  well.  One  platform  may  need  to  be  optimized  for 
hovering  and  agility,  while  another  needs  to  optimized  for  cruising  efficiency.  This  is  where 
classification  of  platforms  comes  into  play.  Secondly,  the  wing  topology  and  wing  kinematics 
must  be  parameterized.  As  seen  in  flapping  wing  species,  wing  topologies  vary  drastically. 
Determining  the  geometric  parameters  that  define  these  various  shapes  is  imperative  for 
successful  optimization.  Not  to  mention,  there  is  the  need  to  parameterize  wing  venation  and 
wing  stiffness.  In  addition,  wing  kinematics  are  more  easily  parameterized  as  having  a  range  of 
plunge,  pitch,  and  lead-lag,  as  well  as  frequencies  and  phase  differences.  Lastly,  all  of  these 
parameters  must  be  constrained  by  a  variety  of  metrics,  whether  it  be  lift/thrust  requirements, 
controllability,  or  manufacturing  capabilities.  Once  all  of  these  are  defined,  trade  studies 
between  these  parameters  can  be  conducted,  and  eventually  various  designs  can  be  optimized 
based  on  a  variety  of  objective  functions  (13-15). 
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5.  Conclusions 


After  using  MBDyn  to  analyze  a  few  flapping  wing  systems,  it  is  clear  that  MBDyn  has  many 
features  that  make  it  a  very  appealing  analysis  tool.  MBDyn’ s  geometry  and  element  definition 
are  extremely  versatile  and  straightforward,  as  is  its  prescribed  dynamic  motion  to  the  system. 
MBDyn  is  capable  of  leveraging  math  libraries  for  eigenvalue  analysis  and  can  even  be  coupled 
with  an  external  CFD  solver.  All  of  these  features  are  important  for  computational  design, 
making  MBDyn  a  nearly  perfect  candidate. 

The  research  did  reveal  some  modeling  issues,  in  particular  the  problem  of  equivalent 
continuum.  Representing  the  wing  surface  as  a  network  of  1-D  beams  is  difficult,  especially 
when  using  an  unstructured  mesh.  Using  polyethylene  material  properties  is  not  necessarily  the 
proper  approach.  As  the  mesh  becomes  finer,  more  beams  are  added  making  the  structure  stiffer. 
In  addition,  to  maintain  the  total  mass  of  the  system  as  the  mesh  becomes  finer,  the  beams  need 
to  become  thinner,  which  again  affects  the  stiffness.  Developing  an  accurate  equivalent 
continuum  method  is  crucial  for  future  modeling  of  FW  systems.  This  problem  does  not  arise 
for  rotor  blade  analysis,  because  airfoil  data  are  used  to  represent  the  blade  surface. 

In  conclusion,  MBDyn  is  a  wonderful  candidate  for  dynamic  analysis  of  aerial  vehicle  platforms 
across  the  board.  Documented  applications  of  using  MBDyn  for  rotor  blade  analysis  prove  its 
capability  in  that  field.  This  research  demonstrated  MBDyn’s  ability  to  model  FW  systems. 
Assuming  the  issue  of  equivalent  continuum  is  addressed,  MBDyn  is  a  strong  candidate  to  be 
used  for  FW  analysis. 
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Appendix  A.  Sample  Mbdyn  Input 


begin:  data; 

problem:  initial  value; 
end:  data! 

begin:  initial  value; 

initial  time:  0. 0000e+00; 

final  time:  1.0000e+00; 

time  step:  1.0000e-03; 

max  iterations:  100; 
tolerance:  1.0000e-05; 

end:  initial  value; 

begin:  control  data; 
structural  nodes: 

+1  #  clamped  node 

+113  #  other  nodes 

rigid  bodies: 

+296  #  mass  of  nodes 

joints : 

+1  #  ground  clamp 

+1  #  total  joint 

beams : 

+296  #  nodal  connectivity 

end:  control  data! 

set:  real  j  =  0. 0000e+00; 

begin:  nodes; 

structural:  1,  static, 

-  1,  0.  ,  0.  , 

eye, 

null, 

null; 

structural:  2,  dynamic, 

5.  9909e-02,  2.  9994e-02, 

eye, 
null, 
null; 

structural. . . 
structural:  114,  dynamic, 


0.  0000e+00, 
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5.  5342e-02,  2.  5161e-02,  0.  OOOOe+OO, 

eye, 
null, 
null; 

end:  nodes; 

set:  real  initial_time  =  0. 0000e+00; 

set:  real  frequency  =  1.2470e+02;  #  radians 

set:  real  amplitude  =  1. 0000e+00; 

set:  real  initial_value  =  0. 0000e+00; 

drive  caller:  100,  sine,  initial_time,  frequency,  amplitude,  forever,  initial_value ; 

reference:  1000, 
null, 
eye, 
null, 
null; 


begin:  elements; 

joint:  10,  clamp,  1,  node,  node; 
joint:  20,  total  joint, 

1, 

position,  reference,  1000,  null, 
position  orientation,  reference, 
rotation  orientation,  reference, 

16, 

position,  reference,  1000,  null, 
position  orientation,  reference, 
rotation  orientation,  reference, 
position  constraint, 

active,  active,  active, 

0. ,  0.  ,  .  001, 
reference,  100, 
orientation  constraint, 

active,  active,  active, 

0.  0,  0.  0,  0.  , 

reference,  100; 


1000,  eye, 
1000,  eye, 


1000,  eye, 
1000,  eye, 


body:  1,  2, 

5.  3988e-06, 

-3.  2671e-04,  -2.  0016e-03,  0.  0000e+00, 

diag,  0. 0000e+00,  7. 4019e-12,  7. 4019e-12; 

beam2:  1, 

2,  null, 

41,  null, 
eye, 

linear  elastic  generic, 

diag,  8.  5322e+02,  8.  5322e+02,  8. 5322e+02,  1.  7674e-04,  7.  1102e-05, 
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1.  0564e-04; 


body:  2,  2, 

4.  2517e-06, 

-2.  0639e-03,  1.  0000e-08,  0.  0000e+00, 

diag,  0. 0000e+00,  6. 0372e-12,  6. 0372e-12; 

beam2 :  2, 

2,  null, 

3,  null, 
eye, 

linear  elastic  generic, 

diag,  1.  5000e+02,  1.  5000e+02,  1. 5000e+02,  2. 5000e-05,  1. 2500e-05, 

1.  2500e-05; 
body.  .  . 


body:  296,  113, 

3.  3209e-06, 

1.2218e-03,  -1.  0517e-03,  0.  0000e+00, 

diag,  0. 0000e+00,  2. 8769e-12,  2. 8769e-12; 

beam2 :  296, 

113,  null, 

114,  null, 
eye, 

linear  elastic  generic, 

diag,  1.  5000e+02,  1.  5000e+02,  1. 5000e+02,  2. 5000e-05,  1. 2500e-05, 

1.  2500e-05; 
end:  elements! 
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Abstract 


Despite  numerous  advances  in  the  field  of  robotics,  robots  remain  astonishingly  poor  at  some 
simple  tasks,  such  as  opening  doors.  In  urban  combat  environments,  doors  are  a  common 
hindrance  to  autonomous  exploration.  As  robots  continue  to  enhance  Soldier  safety  by 
increasing  their  role  in  combat  environments,  autonomous  navigation  of  doors  will  become  ever 
more  important.  Therefore,  the  objective  of  this  project  was  to  implement  and  evaluate  various 
strategies  for  opening  doors  with  a  robotic  manipulator  through  computer  simulation.  A  high- 
fidelity  physics  engine,  Open  Dynamics  Engine,  was  used  to  model  the  Mitsubishi  PA  10 
7-degree-of-freedom  manipulator  in  C++.  To  examine  the  robustness  of  each  strategy  to 
real-world  challenges,  door  characteristics  and  environments  were  varied  in  simulation. 
Environments  analyzed  included  hallways,  rooms,  and  staircase  landings.  Characteristics 
analyzed  encompassed  the  direction  of  the  door’s  motion,  swing  radius,  and  presence  or  absence 
of  auto-closure  devices.  Various  strategies  for  opening  and  closing  doors  were  then 
implemented  in  the  simulation.  These  included  whole -body  motion  using  full  and  selected  joint 
compliance,  real-time  path  planning,  and  equilibrium  point  control.  The  robustness  and  speed  of 
these  strategies  were  then  evaluated  across  the  various  scenarios. 
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1.  Introduction  and  Background 


1.1  Door  Opening  Problem 

In  modem  combat  settings,  thorough  exploration  of  an  environment  is  necessary  to  ensure  safety 
and  mission  success.  Enormous  risks  are  involved  with  the  exploration  of  hostile  areas,  and 
robots  are  frequently  used  to  overcome  these  risks.  Unfortunately,  robots  still  perfonn  poorly 
when  presented  with  tasks  such  as  successfully  navigating  and  manipulating  foreign 
environments. 

The  inability  to  open  doors  speedily  remains  one  of  the  most  common  and  largest  barriers  to  a 
robot’s  usefulness  in  these  situations.  The  primary  reason  for  this  is  the  door’s  constrained 
nature.  Small  errors  in  the  manipulator’s  path  can  create  large  internal  forces,  potentially 
causing  damage  to  the  robot  and/or  doorway  and  resulting  in  a  failure  to  successfully  manipulate 
the  door.  In  addition,  doors  can  be  found  with  widely  varying  features  and  in  widely  varying 
environments.  They  have  different  masses,  swing  radii,  closing  apparatuses,  and  directions  of 
motion.  The  location  of  a  door  may  be  in  an  open  area  or  in  a  confined  space.  Due  to  these 
variations,  a  single  strategy  may  not  be  optimized  for  or  capable  of  handling  the  generic 
problem.  Because  of  this,  different  strategies  may  be  better  suited  to  handle  the  problem  given 
different  scenarios.  Therefore,  this  paper  compares  the  effectiveness  of  several  door  opening 
strategies  in  terms  of  their  speed  and  robustness  to  the  problem  in  various  scenarios. 

1.2  Door  Opening  Strategies 

Compliance  has  been  shown  to  greatly  improve  both  the  speed  and  probability  of  success  of  door 
opening  manipulators  (I).  Allowing  compliance  in  the  system  enables  internal  forces  to  resolve 
themselves.  This,  in  turn,  prevents  binding  and  increases  the  odds  of  success  as  well  as  the 
safety  of  the  action.  Compliance  can  be  achieved  in  a  number  of  ways  including  through 
backdrivability,  the  implementation  of  clutch  mechanisms,  or  inverse  dynamics.  However,  the 
method  of  achieving  compliance  is  outside  the  scope  of  this  paper. 

To  effectively  make  use  of  compliance,  several  strategies  have  been  suggested.  First,  the  door 
might  be  controlled  completely  by  the  mobility  platform,  using  fully  compliant  manipulator 
joints.  This  could  be  implemented  with  a  straight-line  mobility  path  or  with  real-time  course 
correction.  Another  strategy  might  involve  controlling  the  door  with  a  single  manipulator  joint, 
where  the  remaining  joints  are  freely  compliant.  This  would  reduce  the  travel  distance  necessary 
for  the  mobility  platform.  These  strategies  might  also  be  combined  to  handle  scenarios  where 
space  is  limited. 

Finally,  equilibrium  point  control  represents  the  current  state-of-the-art  in  door  manipulation  (2). 
In  this  strategy,  each  joint  is  modeled  as  a  virtual  spring,  enabling  motion  with  low  impedance 
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control.  This  control  method  will  be  implemented  in  simulation  and  will  be  used  for  comparison 
against  the  other  strategies. 

1.3  Door  Opening  Scenarios 

The  optimal  strategy  for  opening  a  door  likely  depends  on  the  scenario  in  which  the  door  is 
placed.  Several  scenarios  were  used  for  testing.  These  included  settings  where  the  door  was 
located  in  a  corner  of  a  room,  at  the  end  of  a  hallway,  and  at  the  top  of  a  landing  for  a  flight  of 
stairs.  Additionally,  we  considered  that  the  most  appropriate  strategy  could  depend  on  the  swing 
radius  of  the  door,  the  mass  of  the  door,  and  the  direction  in  which  the  door  opens.  The  presence 
of  auto-closure  devices  was  also  considered.  This  included  auto-closure  both  as  a  constant  and  as 
a  linearly  increasing  force. 

1.4  Hardware 

The  Mitsubishi  PA10  robot  arm  is  a  commonly  used  7-degree-of- freedom  (7-DOF)  robotic  arm. 
It  contains  seven  rotational  joints,  arranged  as  alternating  pivots  and  hinges.  The  PA  10  attempts 
to  create  a  crude  model  of  the  human  arm  and  is  separated  into  shoulder,  elbow,  and  wrist 
components  (3).  The  PA10  was  selected  for  use  in  this  project  because  of  its  availability,  range 
of  motion,  and  ubiquitous  use  in  robotics  research. 

1.5  Simulation  and  Physics  Engines 

Computer  simulation  is  an  important  part  of  analyzing  potential  solutions  to  robotic  design 
problems.  Simulations  are  a  simple,  cost-effective,  and  safe  alternative  to  field  testing  in 
assessing  the  quality  of  prospective  designs.  To  successfully  model  compliance,  physical 
dependability  is  necessary. 

Physics  engines  are  programming  frameworks  that  use  iterative  methods  to  approximate  rigid 
body  dynamics,  making  high-fidelity  simulation  possible  without  extensive  mathematical 
modeling.  Simulations  that  rely  on  physics  engines  provide  a  high-fidelity  physical  model  while 
using  minimal  computing  power.  In  addition,  performing  experiments  in  simulation  allows 
results  to  be  obtained  much  more  quickly  and  safely  than  conducting  physical  experiments.  For 
these  reasons,  evaluating  various  strategies  and  scenarios  relevant  to  the  door  opening  problem 
through  the  use  of  computer  simulation  is  incredibly  practical. 

1.6  Basis  for  Selecting  Open  Dynamics  Engine 

Open  Dynamics  Engine  (ODE)  was  selected  for  use  based  on  previous  research  comparing 
various  criteria  important  for  use  in  robotic  simulation  through  assessing  both  physical  fidelity 
and  ease  of  development  ( 4 ).  In  the  cited  comparison,  a  weighted  selection  matrix  was  built  that 
contained  a  list  of  important  features  when  using  a  physics  engine.  ODE  was  among  the  physics 
engines  ultimately  selected  for  its  ease  of  use  and  accurate  physics  models.  Based  on  these 
results,  ODE  was  selected  for  use  in  the  door  opening  simulation. 
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1.7  Language  Selection 

ODE  is  implemented  in  the  C  programming  language,  a  highly  popular,  procedural  programming 
language.  Wrappers  for  other  languages  exist,  but  most  of  these  are  poorly  documented,  out  of 
date,  or  incompatible  with  the  current  version  of  ODE.  The  C  implementation  of  ODE  was 
compiled  to  a  dynamically  linked  library  for  use  with  the  project. 

Because  of  the  robust  nature  of  the  door  opening  problem  and  the  breadth  of  scenarios 
constructed,  an  object-oriented  implementation  in  C++  was  used  in  favor  of  a  simple  procedural 
approach.  This  object-oriented  approach  allows  for  the  implementation  of  separate,  independent 
modules,  making  it  easier  to  implement  different  scenarios  without  rewriting  a  significant 
amount  of  code. 


2.  Experiment  and  Calculations 


2.1  Implementation  of  Manipulator 

A  generic  3-degree-of-freedom  (3-DOF)  manipulator  was  built  in  ODE  as  a  proof  of  concept. 
The  manipulator,  shown  in  the  figure  1 ,  contains  three  rotational  degrees  of  freedom,  arranged  as 
a  pivot  joint  and  two  hinge  joints. 


Figure  1.  Generic  3-DOF  manipulator. 

Upon  completion  of  a  generic  manipulator,  the  PA  10  was  coded.  The  PA  10  contains  seven 
rotational  joints,  each  an  alternating  pivot  or  hinge.  The  model,  depicted  in  figure  2,  shows  the 
constructed  PA  10  robotic  arm  approaching  a  generic  door  in  the  comer  of  a  room.  When  the 
PA  10  arm  is  implemented  in  ODE,  it  becomes  possible  to  simulate  and  evaluate  the  different 
scenarios  for  opening  the  door. 
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Figure  2.  Robot  equipped  with  PA10  arm  approaching  a  door. 

2.2  Forward  and  Inverse  Kinematics 

All  the  strategies  used  in  the  door  opening  problem  require  the  simulation  to  attach  the  robotic 
arm  to  a  door  either  before  or  during  the  simulation.  Inverse  kinematics  was  necessary  to 
determine  the  appropriate  joint  angles  needed  to  attach  the  end  effector  to  the  door  (figure  3). 


Figure  3.  PA  10  arm  using  inverse  kinematics  to  reach  a  destination. 
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While  forward  kinematics  seeks  to  determine  the  location  of  an  end  effector  based  on  joint  angle, 
inverse  kinematics  attempts  to  determine  a  set  of  joint  angles  that  a  robotic  arm  needs  to  satisfy  a 
given  end-effector  position  and  orientation. 

Inverse  kinematic  equations  for  the  PA  10  were  evaluated  using  transformation  matrices  of 
Denavit-Hartenberg  parameters.  Generic  inverse  kinematic  equations  from  Robot  Modeling  and 
Control  were  applied  to  the  PA  10  (5).  A  separate  matrix  class  was  constructed  to  supplement  the 
existing  code  in  ODE  and  ensure  that  the  inverse  kinematic  equations  presented  could  be 
successfully  compared  with  the  appropriate  forward  kinematic  solutions. 

In  addition  to  inverse  kinematics,  a  brute  force  method  called  Cyclic  Coordinate  Descent  was 
used  to  navigate  the  end  effector  to  a  desired  position.  This  method  simply  iterates  through  the 
joints  and  attempts  to  move  the  end  effector  closer  to  its  desired  position  by  finding  the  optimal 
angle  for  each  joint.  By  cycling  through  all  the  joints  several  times,  the  arm  typically  becomes  as 
close  as  possible  to  the  target  position.  However,  because  of  its  brute  force  nature,  this  method 
requires  significantly  more  computing  power  than  solving  the  inverse  kinematic  equations  and  is 
less  than  ideal  for  situations  in  which  speed  matters,  such  as  the  door  opening  problem. 

2.3  Strategy  Implementation 

At  the  time  of  this  paper’s  completion,  the  only  strategy  that  had  been  simulated  was  a  free  arm 
controlled  entirely  by  a  mobility  platform  (figure  4).  Given  the  progress  and  modularity  of  the 
code  written,  it  should  not  be  difficult  to  continue  to  implement  various  door  opening  strategies 
into  the  simulation. 


Figure  4.  Using  mobility  platform  to  open  a  door. 
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After  the  complete  construction  and  evaluation  of  the  free  joints  controlled  by  the  mobility 
platform,  a  model  will  be  built  that  includes  individually  controlled  joints.  Data  from  the 
simulated  joint  encoders  will  be  used  to  implement  single-joint  active  control  and  attempt  to 
determine  the  dynamics  of  the  door. 

A  third  scenario  will  use  initial  impulses  rather  than  data  from  joint  encoders  to  attempt  to  model 
the  physics  of  the  door.  The  door  weight  and  friction  will  be  estimated  based  on  these  initial 
impulses. 

Finally,  equilibrium  point  control,  as  described  by  Jain  and  Kemp  (2),  will  be  simulated. 
Equilibrium  control  has  been  shown  to  be  an  effective  method  of  opening  doors  and  drawers  in 
the  past.  This  differs  from  the  previous  scenarios  in  that  the  dynamics  of  the  door  are  not 
considered.  Rather,  “springs”  are  applied  to  arms  to  simulate  the  muscle  found  in  living  beings, 
and  the  links  in  the  arm  are  allowed  to  travel  to  their  equilibrium  points.  These  equilibrium 
points  are  then  modified  to  control  the  ann. 

After  each  scenario  is  implemented,  the  scenarios  will  be  compared  based  on  several  criteria. 
These  criteria  include  the  predicted  speed  at  which  the  scenario  can  be  completed  on  a  physical 
PA  10  manipulator  and  the  force  and  torque  at  the  end  effector.  The  control  schemes  with  online 
dynamics  models  will  be  compared  to  equilibrium  point  control  in  an  attempt  to  determine  the 
benefits  of  using  these  physical  models. 


3.  Results  and  Discussion 


3.1  Physical  Fidelity  of  ODE 

The  physical  fidelity  of  ODE  appears  to  be  great  enough  to  successfully  implement  a  robot, 
manipulator,  and  door  in  a  simulation  environment.  Occasional  joint  instability  was  present  when 
the  robot  was  attached  to  the  ground  with  a  fixed  joint  and  when  the  arm  initially  attached  its 
gripper  to  the  door.  In  both  cases,  the  links  in  the  arm  appeared  to  bounce  beyond  the  constraints 
of  the  joints  connecting  them.  In  addition,  the  door  occasionally  glided  slightly  away  from  its 
axis  when  its  initial  position  and  orientation  were  reset. 

3.2  Results  of  Individual  Strategies 

We  are  continuing  to  work  on  the  results  of  applying  the  various  strategies  discussed  in  section 
2.3.  At  the  time  of  this  publication,  no  results  are  available.  These  results  will  be  published  in  a 
future  document  once  the  project  has  been  completed. 
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4.  Summary  and  Conclusion 


This  paper  discussed  the  difficulties  associated  with  the  door  opening  problem,  several  strategies 
for  opening  a  door  using  a  Mitsubishi  PA  10  robot  arm,  and  the  benefit  of  simulation  and  its 
applicability  to  the  door  opening  problem.  A  PA  10  robotic  arm  was  constructed  in  ODE,  and 
inverse  kinematics  for  the  arm  and  robot  were  implemented.  Several  strategies  for  opening  the 
door  were  tested  and  evaluated,  and  an  online  model  of  the  door’s  dynamics  was  estimated. 
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Abstract 


There  is  growing  interest  to  move  away  from  unguided  munitions  towards  higher-precision, 
guided  munitions.  While  these  smart  projectiles  allow  for  increased  probability  of  hitting  a 
target  and  lower  collateral  damage,  controller  systems  for  these  munitions  are  required  to  be 
small,  rugged,  yet  affordable.  Factors  contributing  to  higher  costs  include  complex  guidance, 
navigation,  and  control  (GNC)  algorithms,  onboard  energetics  (thrusters),  and  large  battery 
power  demands.  One  idea  to  potentially  reduce  these  issues  is  the  concept  of  naturally  roll-stable 
V-tailed  projectiles.  Similar  to  paper  airplanes,  a  roll-stable  projectile  would  be  capable  of 
uprighting  itself  during  flight,  due  to  aerodynamic  loading  in  a  gravity  environment. 
Implementation  of  this  concept  would  provide  numerous  advantages  to  flight  control  systems, 
including  reduced  actuator  burden,  higher  maneuverability,  reduced  sensor  burden,  and 
simplified  GNC  algorithms.  The  scope  of  this  study  is  to  gain  insight  into  the  flight  mechanics 
of  roll-stable  projectiles.  A  flight  dynamic  model  for  asymmetric  V-tailed  projectiles  was 
developed  to  replicate  flight  characteristics  observed  from  initial  testing.  Next,  parametric 
variation  of  V-tail  geometric  parameters  was  performed.  Results  from  these  studies  indicate  that 
naturally  roll-stable  projectiles  can  exist.  Conclusions  of  this  study  are  based  on  aerodynamic 
predictions  and  rigid  6-degree-of-freedom  trajectory  simulation. 
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1.  Introduction 


Projectiles  such  as  mortars  and  artillery  are  used  to  hit  targets;  however,  a  number  of  conditions 
can  cause  rounds  to  miss  an  intended  target.  These  conditions  include  variable  atmospheric 
conditions,  firing  platform  motion,  aiming  errors,  and  manufacturing  inaccuracies  of  the  gun 
tube,  propellant,  and  projectile.  With  the  advent  of  low-cost,  small,  rugged  micro-electro- 
mechanical  systems,  dramatic  reduction  in  dispersion  for  projectiles  equipped  with  a  flight 
control  system  is  possible.  While  these  smart  projectiles  allow  for  increased  probability  of 
hitting  a  target,  as  well  as  lower  collateral  damage,  control  systems  for  these  munitions  are 
required  to  be  small  and  rugged,  yet  affordable. 

Factors  contributing  to  higher  costs  include  complex  guidance,  navigation,  and  control  (GNC) 
components;  onboard  energetics  (thrusters);  and  large  battery  power  demands.  Simplifying  GNC 
systems  effectively  reduces  cost,  but  there  are  many  hurdles  that  result  from  adherence  to 
conventional  projectile  designs.  For  example,  projectiles  are  fired  out  of  rifled  or  smooth-bore 
guns,  which  dictate  the  static  stability.  For  statically  unstable  rounds,  extremely  high  roll  rates 
are  imparted  on  the  projectile  in  a  rifled  gun  barrel.  To  control  these  spin-stabilized  rounds,  the 
GNC  algorithm  is  required  to  process  data  from  high  quality  positional  and  rate  sensors,  and 
coordinate  in  control  mechanism  actuators  with  extremely  fast  response  times.  These  factors 
stress  the  available  technology  and  drive  up  the  cost.  Statically  stable  projectiles  can  mitigate 
these  issues  because  they  ideally  do  not  require  spin  to  fly  efficiently;  however,  sophisticated 
sabot/fin-deployment  designs  are  often  used,  which  again  increase  system  complexity.  Also,  fin 
configurations  are  often  canted  to  induce  a  smaller,  but  still  considerable,  amount  of  spin  so  as  to 
“roll  out”  aerodynamic  and  mass  asymmetries,  thus  reducing  ballistic  dispersion  and  increasing 
precision. 

One  idea  to  potentially  reduce  these  issues  is  to  employ  aerodynamic  and  mass  asymmetries  to 
the  benefit  of  the  control  system.  Specifically,  asymmetric  projectile  airframes  exist  that  do  not 
require  spin.  Similar  in  characteristic  to  the  flight  mechanics  of  a  paper  airplane,  these 
asymmetric  airframes  might  also  have  the  ability  to  naturally  upright  while  in  flight.  If  these 
naturally  uprighting,  roll-stable  projectiles  exist,  and  at  the  same  time  possess  favorable  flight 
characteristics,  then  this  concept  would  provide  numerous  advantages  in  smart  projectile 
applications,  including  reduced  actuator  burden,  higher  maneuverability,  reduced  sensor  burden, 
and  simplified  GNC  algorithms. 

Figure  1  shows  the  concept  of  a  V-tailed  projectile  airframe,  which  is  investigated  in  this  study. 
Similar  to  paper  airplanes,  a  roll-stable  projectile  would  be  capable  of  uprighting  itself  naturally 
during  flight  due  to  aerodynamic  loading  in  a  gravity  environment.  The  scope  of  this  study  is  to 
gain  insight  into  the  flight  mechanics  of  V-tailed  projectiles.  A  flight  dynamic  model  for 
asymmetric  V-tailed  projectiles  was  developed  to  replicate  flight  characteristics  observed  from 
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initial  testing.  Next,  parametric  variation  of  V-tailed  geometric  parameters  was  performed,  and 
finally,  naturally  uprighting  roll-stable  airframes  were  identified. 


Figure  1.  A  3-D  rendering  of  a  V-tailed  airframe  for  proposed  naturally  upright, 
roll-stable  flight  characteristics. 


Figure  2.  Similar  to  paper  airplanes,  V-tailed  projectile  airframes 

utilize  to  flat  aerodynamic  surfaces  located  above  the  center 
of  mass  for  naturally  upright  roll-stable  flight. 
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2.  Projectile  Flight  Dynamic  Model 


2.1  Equations  of  Motion 

The  non-linear  trajectory  simulation  used  in  this  study  is  a  standard  6-degree-of-freedom  (6DOF) 
model  typically  used  in  flight  dynamic  modeling  of  projectiles.  The  6DOF  includes  three 
components  of  the  center  of  gravity  (CG)  position  vector  (x,v,z)  and  three  Euler  projectile 
orientation  angles  (</>,9,y/),  referenced  to  an  inertial  frame.  Refer  to  figure  3  for  an  illustration  of  a 
projectile  possessing  6  degrees  of  freedom. 


Figure  3.  Projectile  axes  (x,y,z)  and  Euler  angles  (<[).0.<|/). 


The  kinematic  and  dynamic  equations  for  the  6DOF  system  are  provided  in  equations  1-4,  as 
derived  by  McCoy  (1 )  and  Carlucci  (2). 
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In  equations  1  and  2,  standard  shorthand  notation  for  trigonometric  functions  is  used:  sa  =  sin(a), 
ca  =  cos(a),  and  ta  =  tan(a).  X,Y,Z  and  L,M,N  are  the  force  and  moment  components  appearing  in 
equations  3  and  4.  Three  translational  velocity  components  ( u,v,w )  and  three  angular  velocity 
components  {p,q,r )  are  related  to  the  time  rate  of  change  of  the  position  vector  (x,y,z)  and 
orientation  angles  y/),  respectively.  Also,  m  and  /  represent  the  projectile  total  mass  and 
inertia  tensor. 


2.2  Gravity  and  Aerodynamic  Modeling 


Projectile  body  forces  were  modeled  using  superposition  with  contributions  from  gravity  (G), 
body  aerodynamics  (A),  and  fins  (F).  Using  superposition,  the  force  components  are 
decomposed,  as  shown  in  equation  5. 
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The  dynamic  equations  are  expressed  in  a  body  fixed  reference  frame,  thus  all  forces  acting  on 
the  body  are  expressed  in  the  projectile  reference  frame.  The  force-acting  due  to  gravity  is 
shown  in  equation  6. 
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The  body  aerodynamic  force-acting  at  the  center  of  pressure  (COP)  of  the  projectile  is  given  by 
equation  7. 
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The  aerodynamic  force  coefficients  are  zero-yaw  drag  ( CXo ),  yaw  drag  (CX2),  and  normal  force 
slope  due  to  angle  of  attack  ( Cna )•  Atmospheric  density  is  expressed  as  p,  and  the  projectile 
diameter  and  total  velocity  are  D  and  V,  respectively.  The  fin  forces  will  be  discussed  in  the  next 
section. 


168 


The  applied  moments  about  the  projectile  mass  center  contains  contributions  from  steady 
aerodynamics  (SA),  unsteady  aerodynamics  (UA),  and  fins  (F). 
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The  moment  components  due  to  steady  aerodynamic  forces  and  control  forces  are  computed  by 
the  cross  product  between  the  distance  vector  from  the  mass  center  to  the  location  of  the  specific 
force  and  the  force,  itself.  The  unsteady  body  aerodynamic  moment  provides  a  damping  source 
for  projectile  angular  motion  and  is  given  by  equation  9. 
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The  aerodynamic  moment  coefficients  are  the  roll  moment  due  to  fin  cant  (Clod),  roll-damping 
moment  (Clp),  and  pitch  damping  moment  ( Cmq )• 


2.3  V-tailed  Fin  Model 


The  fin  forces  are  modeled  as  aerodynamic  lifting  surfaces,  where  the  force  due  to  a  single  fin  is 
modeled  as  a  point  force-acting  at  the  lifting  surface  aerodynamic  center  of  pressure.  The 
moment  generated  due  to  lift  about  the  fin  quarter-chord  is  neglected  in  this  study.  Fin 
orientation  on  the  projectile  is  defined  as  a  set  of  three  body-fixed  rotations  about  the  projectile 
reference  frame.  Starting  with  the  canard  axis  aligned  with  the  projectile  body  axis,  the  ith  fin  is 
rotated  about  the  is  axis  by  the  azimuthal  angle  (</>,),  then  about  the  resulting  intermediate  /c-axis 
by  the  sweep  angle  (yl j.  A  third  angle,  the  fin  cant  angle  (Sl),  will  be  addressed  when  defining 
local  fin  angle  of  attack.  If  all  three  fin  angles  are  zero,  the  lifting  surface  is  in  the  fg  -  Jb  plane. 
The  transformation  from  the  ith  fin  reference  frame  to  the  parent  projectile  body  axis  is  given  by 
equation  10. 
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Using  equation  10,  the  body  referenced  air  velocities  u,  v,  and  w  of  equation  3  can  be 
transformed  into  the  ith  fin  reference  frame  as  w„  v„  and  vv;-,  to  obtain  local  angle-of-attack  for 
aero  calculations.  The  vector  components  rxi,  ryi,  and  rzi  represent  the  distances  between  the 
projectile  CG  and  the  ith  fin. 
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From  this,  the  fin  forces  can  be  modeled  as  a  function  of  the  local  angle-of-attack,  which  for  the 
ith  lifting  surface  is  given  by  equation  12. 


or,  =  8j  +  tan 
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The  total  aerodynamic  force  components  generated  by  the  two  V-tailed  fins  and  expressed  in  the 
projectile  body  frame  is  given  by  equation  13. 
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The  moment  components  generated  by  the  fins  on  the  projectile  body  are  computed  by  the  cross 
product  between  the  distance  vector  from  the  CG  to  the  location  of  the  ith  fin  force,  previously 
described.  Then  components  for  both  fins  are  summed  to  obtain  the  total  moment. 


The  aerodynamic  coefficients  Cu  and  Cm  are  modeled  according  to  the  following  expressions. 


Cl.  —  C,Aal  +  Cf3al  +  Cl5al 

(14) 

Cn  =  CD0  +  CD2a~  +  C,C,2 

(15) 

The  ith  fin  COP  location  and  the  aerodynamic  coefficients,  Cu  and  CDi,  are  again  predicted  as 
functions  of  Mach  number  and  surface  angle-of-attack.  The  fin  transverse  COP  location  was 
assumed  to  act  in  the  transverse  CG  location  of  the  fin. 


Along  with  Mach  number  and  angle-of-attack,  the  aerodynamic  fin  coefficients  of  equations  14 
and  15  are  also  dependent  upon  fin  geometry.  Six  parameters  were  identified  to  be  important  in 
evaluating  V-tailed  airframes  for  performance:  fin  width  (Wfm),  fin  length  (Lfm),  sweep  angle 
(SA),  V-tailed  angle  (0V),  a  fin  offset  parameter  (h),  and  fin  cant  angle.  An  illustration  of  these 
parameters  is  shown  in  figure  4. 
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Side  View: 


Rear  View: 


l -fin 


Figure  4.  Illustration  of  V-tailed  fin  geometry  parameter  space. 

2.4  Simulation  Implementation 

The  dynamic  equations  given  by  equations  1-4  are  numerically  integrated  forward  in  time  using 
a  4th-order,  fixed-step  Runge-Kutta  algorithm.  The  mass,  CG  location,  and  inertia  tensor  are  all 
assumed  to  be  constant  throughout  the  duration  of  the  flight.  The  COP  location  and  all 
aerodynamic  coefficients  in  equations  7,  9,  14,  and  15  are  estimated  using  aero-prediction 
techniques  as  functions  of  Mach  number  and  angle-of-attack.  During  simulation,  these 
coefficients  are  linearly  interpolated. 


3.  Results  and  Discussion 


3.1  Projectile  Solid  Modeling 

Customized  projectiles  were  modeled  in  SolidWorks  to  obtain  exterior  geometries  and  physical 
properties.  These  data  were  used  to  estimate  aerodynamic  coefficients.  Note  that  for  all  6DOF 
simulations,  the  mass,  inertia,  and  CG  location  were  held  constant  for  varied  fin  configurations. 
This  is  because  even  though  modifying  fin  parameters  will  consequentially  change  physical 
parameters,  the  changes  are  expected  to  be  reasonably  small,  since  the  fin  mass  and  volume  are 
much  smaller  than  the  projectile  body. 

3.2  Initial  Flight  Experiments 

Initial  flight  experiments  were  conducted  on  asymmetric  V-tailed  projectiles  with  different  fin 
configurations.  A  low-speed  airgun  was  used  as  a  test  rig,  where  roll  performance  was  evaluated 
by  line-of-site  observation.  Projectile  airframes  were  considered  naturally  roll-stable  if  the  roll 
oscillation  damped  throughout  flight  when  perturbed  by  an  external  influence  (tip-off  conditions, 
wind  gusts,  etc.).  From  these  experiments,  no  airframes  were  seen  to  be  roll-stable,  but  otherwise 
had  good  flight  characteristics — meaning  small  angles-of-attack  were  experienced  through  flight. 

Additional  tests  were  performed  with  video  cameras  placed  in  airframes  and  pusher  test  articles 
to  observe  the  roll  motion  of  the  projectile.  The  results  of  these  experiments  show  that  the 
projectile  rolls  through  multiple  cycles  with  a  roll-dependent  angular  velocity  that  gradually 
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increases  in  frequency  through  flight.  Based  on  video  footage,  projectiles  were  observed  to 
slowly  roll  counter-clockwise  and  then,  when  viewed  from  behind,  change  direction  as  the 
projectiles  became  roll-unstable.  Also,  projectile  states  at  muzzle  exit  were  estimated  to  be  used 
in  6DOF  simulations,  as  shown  in  table  1 . 


Table  1.  Projectile  initial  conditions  to  be  used  in  a  nominal  6DOF  simulation. 
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3.3  6DOF  Model  Validation 

One  of  the  projectile  test  articles  was  chosen  to  be  a  baseline  fin  geometry  for  purposes  of  6DOF 
simulation:  Lfm= 2.42847  (cal),  Wfm=0. 709695  (cal),  &4=54.0  (deg),  h= 0.00  (cal),  <9V=100  (deg), 
and  <5—1.0  (deg),  where  1.0  (cal)  is  equivalent  to  the  characteristic  projectile  diameter.  Airgun 
tip-off  was  estimated  (i.e.,  initial  pQ,  q0,  and  rQ  states)  until  the  roll  characteristics  were  observed. 
The  plots  shown  in  figure  5  compare  single  trajectory  results  for  a  projectile  with  a  tip-off  of  p„= 
0.03  (rad/s),  q0=  -0.22  (rad/s),  and  rQ=  0.052  (rad/s),  and  without  tip-off. 
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Figure  5.  Selected  6D0F  simulation  plots  of  the  baseline  V-tailed  projectile  with  and  without  estimated  tip- 
off  effects  including  (a)  altitude  vs.  time,  (b)  deflection  vs.  time,  (c)  roll  rate  vs.  time,  (d)  roll  angle 
vs.  time,  (e)  pitch  angle  vs.  time,  (f)  yaw  angle  vs.  time,  (g)  pitch  aerodynamic  angle-of-attack  vs. 
time,  (h)  aerodynamic  angle  of  sideslip  vs.  time,  and  (i)  pitch  angle  vs.  angle  of  sideslip. 

In  figures  5(a)  and  5(b),  trajectory  results  show  that  the  projectile  with  tip-off  has  yaw 
characteristics.  Figures  5(d)  -5(e)  show  Euler  angles,  where  figure  5(c)  is  the  roll  rate  time 
history  of  the  projectile  throughout  flight.  As  can  be  seen  from  both  the  roll  angle  and  roll  rate, 
the  projectile  initially  begins  to  slowly  roll  in  the  positive  direction  but  then  changes  direction 
and  starts  to  roll  in  a  cyclic  fashion.  The  roll  rate  increases  to  a  nearly  constant  frequency  until 
late  in  the  flight  when  oscillations  are  observed.  Note  that  after  approximately  2.5  s  of  flight,  the 
oscillatory  behavior  is  observed  in  all  three  Euler  angles.  Figures  5(g)  and  5(h)  show  projectile 
aerodynamic  angles’  time  history.  Initial  oscillations  are  observed  in  the  pitch  angle  but  not  in 
the  sideslip  angle.  As  time  progresses,  these  angles  start  to  develop  growing  oscillatory 
characteristics,  which  represents  decreased  flight  performance. 
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3.4  Mechanics  Driving  Fundamental  Projectile  Behavior 

Several  simulations  were  computed  to  study  the  effects  of  gravity  and  fin  aerodynamics.  The 
contributing  forces  due  to  fin  geometry  are  desired  to  keep  the  projectile  roll-stable  in  flight, 
where  gravity  uprights  the  projectile  during  flight.  Throughout  the  trajectory,  gravity  causes  the 
projectile  to  follow  a  near  parabolic  path  in  the  downrange  altitude  plane.  From  this,  a  small  lag 
exists  between  the  spin-axis  and  the  CG  motion,  inducing  projectile  angle-of-attack.  The  roll- 
instability  is  produced  from  the  coupling  of  this  gravity-induced  small  angle-of-attack  and  the 
fins.  Inspection  of  equations  11-13  shows  that  gravity-induced  q  and  r  terms  produce 
asymmetric  fin  forces  on  opposite  sides  of  the  V-tail,  which  can  drive  the  roll-instability. 

Figure  6  shows  selected  results  from  a  study  investigating  the  effects  of  gravity  and  fin  aero  on  a 
projectile  body.  Because  this  study  is  purely  academic,  no  tip-off  effects  were  modeled.  Also  in 
the  absence  of  fins,  the  CG  was  artificially  moved  forward  to  eliminate  static  instabilities 
commonly  observed  in  finless  projectiles.  This  does  not  affect  motion  in  the  roll  direction. 
Shown  in  Figures  6(a)  and  6(b),  the  V-tailed  projectile  in  the  presence  of  a  gravity  environment 
becomes  roll-unstable,  even  in  the  absence  of  tip-off;  however,  the  finless  projectile  cases  do  not 
exhibit  any  motion  in  the  roll  direction.  Figure  6(c)  shows  that  the  introduction  of  gravity  to  the 
finless  projectile  introduces  oscillations  in  the  vertical  plane,  and  when  fins  are  present,  a  small 
offset  pitch  angle  is  induced.  Once  the  finned  projectile  reaches  the  trajectory  apogee,  the 
instability  is  observed  (after  approximately  4.3  s  of  flight).  Figure  6(d)  shows  that  the  instability 
in  the  V-tailed  projectile  is  not  isolated  to  the  vertical  pitch  plane,  but  is  coupled  in  the  lateral 
yaw  plane.  This  confirms  the  predictions  of  equations  1 1-13. 
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Figure  6.  Selected  6D0F  simulation  plots  showing  how  of  gravity  and  fin  aero  effects  on 
projectile  roll  stability:  (a)  roll  angle  vs.  time,  (b)  roll  rate  vs.  time,  (c)  pitch 
angle  vs.  time,  and  (d)  yaw  angle  vs.  time. 

3.5  V-tailed  Fin  Parametric  Study 

To  identify  possible  roll-stable  airframes,  fin  parameters  were  parametrically  varied  from  the 
baseline  fin  geometry  to  map  out  roll-stability  performance.  These  parameters  include  Wfm,  Lfm, 
h,  6V,  and  d.  For  simplicity,  the  fin  sweep  angle  was  allowed  to  vary  freely  with  changes  in  Ljin 
and  Wfm,  which  is  a  valid  assumption  for  the  low  flight  speeds  observed  in  simulation.  New 
estimated  aero  predictions  were  computed  based  off  of  new  fin  geometry  prior  to  each 
parametric  study.  The  total  projectile  mass  and  inertial  properties  were  not  recalculated, 
however.  Observations  of  these  parametric  studies  are  summarized  below  in  table  2. 
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Table  2.  Results  and  summary  of  fin  parametric  studies. 


Fin  Parameter 

Trends  and  Observations 

// 

■  Without  tip-off  effects,  increasing  the  fin  offset  parameter  from  0.0  (cal)  to  0.5 
(cal)  (the  top  of  the  projectile)  can  delay  roll-instabilities  in  flight  but  not 
eliminate  them. 

■  If  tip-off  effects  are  considered,  this  parameter  does  not  dampen  roll-instability 
or  promote  projectile  uprighting. 

Wfin 

■  Increasing  fin  width  from  the  baseline  was  not  observed  to  positively  affect 
projectile  performance. 

Lf,„ 

■  Increasing  fin  length  from  the  baseline  was  not  observed  to  positively  projectile 
performance.  Degradation  due  to  variation  of  Lfm  is  not  nearly  as  sensitive  as 

wfln. 

0 

■  Decreasing  9V  increases  stability.  For  60.0  (deg)  >  9V  >  55.0  (deg),  roll-stable, 
naturally  uprighting  airframes  were  observed. 

■  Increasing  h  and  decreasing  S  both  increased  the  0V  roll-stable  window  form 

55.0  (deg)  up  to  approximately  75.0  (deg). 

8 

■  There  exists  a  range  (from  -1.0  (deg)  to  -3.5  (deg))  for  which  the  fin  cant  angle 
can  be  adjusted  to  make  the  V-tail  appear  to  be  roll-stable  throughout  flight  but 
no  naturally  uprighting  airframes  were  observed,  especially  in  the  presence  of 
tip-off  effects. 

■  As  the  fin  can  angle  increases,  the  steady  state  aerodynamic  angle-of-attack 
increases  twice  as  much. 

3.6  Roll-Stable  Projectile  Flight  Characteristics 

Given  the  knowledge  obtained  from  the  parametric  studies,  several  different  fin  geometries  were 
simulated  to  determine  roll-stable  airframes.  As  previously  discussed,  airframes  were 
determined  to  be  naturally  uprighting  and  roll-stable  if,  for  any  set  of  tip-off  conditions,  the  Euler 
roll  angle  was  observed  not  to  fall  into  the  common  cyclic  rolling  motion  but,  instead,  oscillate  in 
a  dampening  fashion  about  the  upright  orientation  of  (j)  =  0.0  (deg).  Additionally,  the  robustness 
of  a  roll-stable  airframe  is  assessed  by  simulating  projectile  fly-outs  for  varying  tip-off  intensities 
and  directions. 

Figure  7  shows  results  of  a  tip-off  analysis  for  a  roll-stable  configuration.  The  trajectory  results 
in  figures  7(a)-7(b)  show  that  tip-off  effects  influence  lateral  displacement  of  the  projectile,  but 
that  altitude  trajectory  is  virtually  unchanged.  Figures  7(d)  -7(e)  show  Euler  angles,  where 
figure  7(b)  is  the  projectile  roll  rate.  The  roll  rate  is  observed  to  mitigate  tip-off  effects  very 
early  into  flight  and  remains  very  low  through  the  remainder  of  flight.  As  a  result,  the  roll  angle 
oscillations  are  reasonably  small  and  centered  about  (f>  =  0.0  (deg)  -  the  upright  position.  Tip-off 
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variation  does  appear  to  affect  the  yaw  angle  in  figure  7(e)  but  not  the  pitch  angle  in  figure  7(f). 
Figures  7(g)-7(h)  show  that  all  projectile  aerodynamic  angles  dampen  throughout  flight,  which 
is  favorable.  Since  gravity  is  a  known  source  of  potential  energy,  the  roll  rate  and  angle  plots 
indicate  that  the  fin  aero  has  a  damping  effect  in  the  roll  plane  for  certain  fin  configurations. 
Thus,  naturally  uprighting,  roll-stable  projectile  airframes  can  exist. 


Figure  7.  Selected  6DOF  simulation  plots  of  a  roll-stable  V-tailed  projectile  for  varied  tip-off  conditions 

including  including  (a)  altitude  vs.  time,  (b)  deflection  vs.  time,  (c)  roll  rate  vs.  time,  (d)  roll  angle 
vs.  time,  (e)  pitch  angle  vs.  time,  (f)  yaw  angle  vs.  time,  (g)  pitch  aerodynamic  angle-of-attack  vs. 
time,  (h)  aerodynamic  angle  of  sideslip  vs.  time,  and  (i)  pitch  angle  vs.  angle  of  sideslip. 


4.  Conclusion 


This  research  indicates  that  naturally  uprighting,  roll-stable  airframes  can  exist  with  favorable 
flight  characteristics  and  should  be  considered  further  for  possible  implementation  in  precision- 
guided  applications.  Roll- instabilities  observed  in  flight  experiments  for  a  chosen  baseline 
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airframe  were  able  to  be  computationally  reproduced,  where  fins  were  modeled  as  flat  surfaces. 
Next,  fin  parameters  were  varied  from  the  baseline  airframe  to  determine  how  each  affects 
projectile  performance.  Increasing  the  h  positively  affected  the  roll-stability  positively,  but  with 
limited  influence  when  tip-off  was  introduced.  Increasing  both  Wfm  and  Lfin  were  not  observed  to 
encourage  roll  stability;  however,  0V  and  5  were  both  observed  to  significantly  promote  roll- 
stability  within  certain  limited  regions.  Using  these  trends,  several  roll-stable  airframes  were 
able  to  be  simulated.  Further,  some  roll-stable  airframes  were  observed  to  be  fairly  robust  to 
variation  in  tip-off  effects.  Conclusions  of  this  study  are  based  on  aero-prediction  techniques  and 
rigid  6DOF  trajectory  simulation. 


5.  Future  Work 


A  second  round  of  flight  experiments  is  currently  being  performed  to  demonstrate  and  confirm 
simulations  that  naturally  upright  roll-stable  airframes  exist.  Upon  successful  flight  testing,  there 
are  several  directions  that  could  be  explored.  For  example,  a  more  sophisticated  test  rig  could  be 
developed  for  more  useful  results  (i.e.,  better  6DOF  validation,  etc.).  Also,  6DOF  simulations 
using  control  mechanisms  could  be  developed  to  determine  projectile  controllability  and 
corresponding  flight  characteristics  while  under  control.  Further  attention  should  also  be 
invested  into  obtaining  better  aero-predictions  using  coupled  CFD-6DOF  methods  proposed  by 
Sahu  ( 6 ),  so  that  complicated  flow  effects  such  as  fin  shadowing  and  flow  separation  can  be 
accounted. 
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Abstract 


Integration  of  synthetic  and  biomaterials  into  novel  hybrid  systems  has  potential  utility  in  a 
variety  of  Army  applications.  Chitosan  is  a  biocompatible  and  biodegradable  material  that  can 
be  produced  in  large  quantities.  This  project  focuses  on  exploring  alternative  coupling  methods 
and  chemistries  for  chitosan  functionalization.  As  a  proof-of-principle,  chitosan  is  being  grafted 
with  polyethylene  glycol  (PEG),  polypropylene  oxide  (PPO),  or  dodecane  to  alter  the  pure 
chitosan  wettability.  Preliminary  Fourier  transform  infrared  spectroscopy  (FTIR)  results  and 
contact  angle  measurements  provide  contradictory  data  regarding  the  degree  of  functionalization. 
This  contradiction  is  currently  being  explored  through  repeated  experimentation  and  may  be  due 
to  the  preferential  segregation  of  low  surface  tension  groups  to  the  air  interface  during  film 
casting.  Angle-resolved  x-ray  photoelectron  spectroscopy  (XPS)  and  synchrotron  near  edge 
x-ray  absorption  fine  structure  (NEXAFS)  will  be  explored  to  examine  the  surface  chemistries. 
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1.  Introduction/Background 


The  incorporation  of  controlled  functionality  of  synthetic  chemistries  into  biological  systems  has 
the  potential  to  impact  a  variety  of  Army-relevant  applications  including  tissue  scaffolding, 
chemical  or  biological  decontamination,  battlefield  wound  treatment,  and  environmental 
remediation.  A  key  aspect  of  many  of  these  applications  is  the  ability  to  tailor  the  interaction  of 
the  hybrid  material  to  its  intended  environment.  For  example,  the  synthetic  material  could  illicit 
a  response  to  a  specific  analyte,  such  as  for  controlled  release  of  a  therapeutic  agent  or  chemical 
countermeasure.  Conversely,  the  tailored  interaction  could  selectively  sequester  a  specific 
species  to  aid  decontamination.  In  most  of  these  applications,  a  biologically  compatible  material 
is  desirable  so  that  there  is  minimal  impact  on  the  biological  system. 

A  promising  biocompatible  scaffold  material  is  chitosan.  Chitosan  is  a  deacetylated  form  of 
chitin,  which  is  a  material  found  in  the  shell  of  many  crustaceans.  Chitosan  is  one  of  the  most 
abundant  biodegradable  materials,  making  it  an  excellent  candidate  for  use  in  sensitive  biological 
environments.  Additionally,  its  abundance  potentially  allows  any  developed  technology  to  scale- 
up  to  useful  quantities  while  maintaining  a  low  cost.  Chitosan  is  a  linear  polysaccharide 
composed  of  D-glucosamine  and  N-acetyl-D-glucosamine  (figure  1).  Due  to  its 
polyglucosamine  structure,  which  is  a  common  dietary  fiber  similar  to  cellulose,  chitosan  is 
generally  considered  safe  for  human  consumption  (7). 


°„ 
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Figure  1.  Chemical  structure  of  chitosan. 

Chitosan’s  chemical  structure  contains  both  hydroxyl  and  amine  functionalities  that  can  be  used 
for  modification  with  synthetic  molecules.  The  focus  of  this  project  is  to  change  the 
hydrophobicity/hydrophilicity  of  chitosan,  primarily  monitored  by  contact  angle.  When  a  water 
droplet  is  placed  onto  a  substrate,  it  forms  a  bead;  the  bead’s  shape  is  dependent  on  the 
interaction  between  the  water  and  the  surface.  By  modifying  the  surface  with  a  hydrophilic 
material,  the  surface  becomes  attracted  to  water  and  decreases  the  contact  angle.  In  contrast,  a 
hydrophobic  modification  would  make  the  newly  functionalized  surface  repel  water  and  the 
contact  angle  would  increase.  This  effect  is  most  commonly  observed  during  the  process  of 
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waxing  a  vehicle;  the  addition  of  the  hydrophobic  wax  causes  the  water  to  bead  up  and  roll  off 
the  car.  This  project  explores  a  range  of  possible  combinations  of  hydrophobic  and  hydrophilic 
polymer  attachments  to  chitosan  as  a  means  of  controlling  the  wettability,  adhesion,  and 
solubility  in  various  media. 

Three  polymer  materials  were  evaluated  for  their  impact  on  wettability  when  coupled  with 
chitosan.  Polyethylene  glycol  (PEG)  is  a  hydrophilic  polymer  that  is  commonly  used  in 
biological  environments  including  therapeutics,  personal  hygiene  products,  and  biological 
implant  coatings.  It  is  soluble  in  water  partially  due  to  its  ability  to  form  hydrogen  bonds  (figure 
2a).  Polypropylene  oxide  (PPO)  is  very  similar  to  PEG  except  it  has  an  additional  methyl  group 
attached  to  the  molecule  (figure  2b).  This  subtle  change  to  the  polymer  backbone  is  enough  to 
disrupt  the  hydrogen  bonding  character  of  the  polymer  and  dramatically  change  its  water 
solubility.  Dodecane  consists  of  a  long  carbon  chain  that  exhibits  a  high  degree  of 
hydrophobicity  due  to  the  large  number  of  methylene  groups  present  (figure  2c). 


Figure  2.  Chemical  structures  of  (a)  polyethylene  glycol,  (b)  polypropylene  oxide  and  (c)  dodecane. 


To  alter  the  wettability  and  adhesion  of  chitosan,  the  hydrophobic  and  hydrophilic  materials  must 
be  functionalized  onto  chitosan  so  that  the  chemistry  of  hybrid  material  changes.  A  major 
drawback  of  chitosan  is  its  lack  of  solubility  in  most  solvents,  including  water,  making  it 
extremely  difficult  to  modify  chitosan’s  chemistry.  To  obtain  an  aqueous  solution  of  chitosan, 


the  water  must  be  acidic,  limiting  the  available  chemistries  for  coupling  (2).  As  a  result,  the 


chitosan  has  been  functionalized  to  enhance  its  solubility  in  various  solvents.  Chitosan 
functionalization  has  been  attempted  through  a  variety  of  synthetic  approaches  (2).  These 
approaches  typically  require  complex  synthetic  methodologies  to  obtain  a  large  fraction  of 
coupling  on  the  chitosan  backbone.  For  this  project,  it  may  not  be  necessary  to  obtain  a  high 
degree  of  coupling  if  a  critical  degree  of  functionalization  can  be  identified  that  produces  the 
required  property  changes. 

The  main  focus  of  this  project  is  to  determine  the  viability  of  alternative  reaction  conditions  for 
the  functionalization  of  chitosan  to  ultimately  alter  its  wettability.  The  amine  group  on  the 
chitosan  is  a  candidate  for  a  simple  Michael  addition  by  reacting  with  a  carbon-carbon  double 
bond.  The  Michael  addition  is  facilitated  by  a  methyl  methacrylate  terminated-PEG  (figure  3). 
The  amine  group  can  also  react  with  an  epoxy  functionalized  polymer.  The  epoxy-amine  is  a 
robust  reaction  found  in  many  commercially  available  adhesives.  This  reaction  is  initially 
investigated  using  a  mono-functional  epoxy  dodecane,  but  can  be  extended  to  epoxy- 
functionalized  PEG  if  the  results  are  promising  (figure  3). 
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For  many  of  the  potential  applications,  a  solution  of  dispersed  chitosan  will  hinder  the  ability  to 
effectively  deploy  or  contain  its  intended  target.  As  a  result,  the  project  also  investigates  the 
ability  to  use  the  previously  described  chemistries  to  attach  the  modified  chitosan  to  scaffolding 
materials  with  useful  geometries.  For  example,  the  epoxy  chemistry  has  been  implemented  into 
the  production  of  foam  materials  that  have  a  large  surface  area.  In  addition,  the  modified 
chitosan  can  be  placed  onto  a  monolithic  material  to  enhance  the  barrier  properties.  It  is 
anticipated  that  these  approaches  can  be  readily  scaled-up  to  accommodate  current  U.S.  Army 
needs. 


2.  Experimental  Setup 


2.1  Chitosan  Functionalized  with  PEG 

Using  similar  protocols  found  in  the  literature,  2%  chitosan  solutions  were  dissolved  in  acidified 
water  using  0.2%  acetic  acid.  PEG  was  then  added  into  the  solution;  the  amount  of  PEG  was  one 
of  the  following:  25%,  50%,  or  75%  by  weight  of  chitosan  in  the  solution.  Blends  of  the  two 
components  were  created  by  mixing  the  solutions  for  an  hour  and  then  directly  casting  the 
solutions  into  silicon  molds.  Reacted  solutions  were  also  made  by  constantly  stirring  the  mixture 
at  60  °C  and  500  rpm  for  4  h  under  nitrogen  atmosphere  (5).  Alternatively,  1%  chitosan/2.0% 
acetic  acid  solutions  with  the  same  amounts  of  PEG  were  produced.  The  solutions  were 
constantly  stirred  at  60  °C  and  500  rpm  for  4  h  under  nitrogen  atmosphere. 

All  PEG  samples  were  washed  with  the  following  method,  except  the  2%  chitosan  PEG  blends. 
First,  10%  w/w  sodium  hydroxide  (NaOH)  was  used  to  precipitate  the  polymer  out  of  the 
solutions.  The  obtained  polymer  was  submerged  in  methanol  and  centrifuged  to  remove  any 
unreacted  molecules  (3).  This  washing  protocol  was  repeated  three  times,  producing  a 
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condensed  pellet.  The  pellet  was  then  redissolved  into  acidified  water  and  cast  into  silicone 
molds  to  form  thin  films  used  for  contact  angle  measurements  and  Fourier  transform  infrared 
spectroscopy  (FTIR). 

2.2  Chitosan  Functionalized  with  PPO 

Various  amounts  of  PPO  (25%,  50%,  and  75%  by  weight  of  chitosan)  were  stirred  into  2% 
chitosan/0.2%  acetic  acid  solutions  for  1  h  to  create  2%  chitosan/PPO  blends.  Additionally,  1% 
chitosan/2.0%  acetic  acid  solutions  were  also  mixed  with  25%,  50%,  and  75%  PPO  at  1000  rpm 
and  room  temperature  for  approximately  24  h. 

2.3  Chitosan  Functionalized  with  Epoxy  Dodecane 

At  this  point,  1%  chitosan  solutions  were  dissolved  in  30%  acidified  water,  using  0.2%  acetic 
acid,  and  70%  ethanol.  Then,  50%  epoxy  dodecane  (half  of  the  weight  of  chitosan)  was  added 
into  the  chitosan  solution.  The  solutions  were  constantly  stirred  at  90  °C  and  500  rpm  for  4  h 
under  nitrogen  atmosphere.  The  same  NaOH  washing  protocol  was  used  to  extract  the  polymer 
and  cast  it  into  films. 


3.  Results/Discussion 


3.1  Pure  Chitosan 

Contact  angle  measurements  were  used  to  quantify  the  wettability  of  the  chitosan-functionalized 
films.  The  instrument  places  a  precise  droplet  of  water  onto  the  substrate,  i.e.,  chitosan,  and  is 
equipped  with  a  camera  to  capture  the  droplet  and  measure  the  contact  angle.  Pure  chitosan 
exhibits  a  contact  angle  of  1 10°  as  seen  in  figure  4. 


Figure  4.  Contact  angle  of  pure  chitosan  (110°). 

FTIR  was  used  to  monitor  the  extent  of  chitosan  with  the  desired  coupling  agent.  FTIR 
measures  the  absorbance  of  the  material  at  a  range  of  wavelengths  that  can  be  assigned  to 
specific  functional  groups.  Figure  5  contains  a  typical  FTIR  scan  of  a  pure  chitosan  film.  The 
key  peaks  include  the  amine  functionality  at  1600  cm  1  and  the  carbon-hydrogen  (carbogen) 
bond  stretching  around  2900  cm-1. 


190 


Figure  5.  FTIR  scan  for  pure  chitosan  with  carbon-hydrogen  bond  stretching  occurring  around  2900  cm 
and  amine  groups  around  1600  cnT1. 

3.2  PEG  Samples 

3.2.1  2%  Chitosan/0.2%  Acetic  Acid  Samples 

To  determine  the  approximate  amount  of  PEG  required  to  alter  the  hydrophilicity  of  chitosan, 
films  consisting  of  chitosan-PEG  blends  were  produced  where  the  two  components  were  mixed 
but  not  reacted.  The  2%  chitosan  blends  revealed  a  significant  change  in  contact  angle  with  the 
addition  of  PEG.  Samples  containing  25%  PEG  and  50%  PEG  samples  exhibited  significant 
increases  in  hydrophilicity  with  contact  angles  of  39.33°  ±  3.93°  (Figure  6a)  and  31.96°  ±  9.34° 
(figure  6b),  respectively.  In  contrast,  the  75%  PEG  blend  became  more  hydrophobic,  with  a 
contact  angle  of  87.79°  ±  15.74°  (figure  6c). 


Figure  6.  Contact  angle  of  2%  chitosan/0.2%  acetic  acid  and  a  (a)  25%  PEG  blend  (39°),  (b)  50%  PEG 
blend  (32°),  and  (c)  75%  PEG  blend  (88°). 
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The  FTIR  scans  of  the  three  blends  support  the  presence  of  PEG  in  the  samples  (figure  7).  The 
PEG  carbon-oxygen  double  bond  (carbonyl)  peak  is  present  at  around  1710  cm-1.  In  addition, 
the  peak  around  2900  cm  1  increases  due  to  the  additional  hydrocarbons  in  the  PEG  backbone. 
The  peak  areas  are  also  consistent  with  the  PEG  concentration  in  the  blends  where  an  increasing 
PEG  concentration  leads  to  an  increased  peak  area.  The  increased  PEG  concentration  in  the  films 
produces  a  contact  angle  decrease  when  the  concentration  of  PEG  in  the  blend  increases  from 
25%  PEG  (39°)  to  50%  PEG  (32°).  However,  the  FTIR  does  not  indicate  why  the  blend 
containing  75%  PEG  would  exhibit  a  decreased  hydrophilicity.  It  is  possible  that  the  high  PEG 
concentration  leads  to  phase  separation  and  inconsistencies  in  the  film  with  regions  of  high  and 
low  PEG  content,  which  may  explain  the  large  standard  deviation.  It  is  also  possible  that  the 
higher  concentration  leads  to  a  preferential  segregation  of  low  surface  tension  groups  to  the  air 
interface  during  film  casting.  To  fully  explore  these  possibilities,  angle-resolved  x-ray 
photoelectron  spectroscopy  (XPS)  or  synchrotron  near  edge  x-ray  absorption  fine  structure 
(NEXAFS)  could  be  used  to  examine  the  surface  chemistries. 


Figure  7.  FTIR  scans  of  the  2%  chitosan/0.2%  acetic  acid  and  25%  PEG  (purple),  50% 

PEG  (green),  and  75%  PEG  (red). 

To  covalently  attach  the  PEG,  the  sample  mixtures  were  heated  to  60  °C  to  react  for  4  h.  The 
samples  were  thoroughly  washed  to  remove  any  unbound  PEG  prior  to  FTIR  analysis.  As  a 
result,  the  presence  of  carbonyl  and  increase  of  the  carbon-hydrogen  peaks  indicates  that  the 
PEG  is  covalently  bonded  or  strongly  attached  to  the  chitosan. 

In  the  2%  chitosan/0.2%  acetic  acid  mixtures  containing  25%  and  75%  PEG,  no  carbonyl  peak 
around  1710  cm'1  was  present  (figure  8a  and  8b,  respectively).  This  indicates  that  the  Michael’s 
addition  reaction  did  not  occur,  and  that  after  the  washing,  all  that  remained  was  pure  chitosan. 
Additionally,  both  samples  had  approximately  the  same  contact  angle  results  as  pure  chitosan 
with  angles  of  109°  (figures  9)  and  104°,  respectively,  indicating  that  the  reaction  did  not  occur. 
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Figure  8.  FTIR  scans  of  2%  chitosan/0.2%  acetic  acid  and  (a)  25%  PEG  and  (b)  50%  PEG 
illustrating  the  similarity  between  the  25%  and  50%  samples  and  pure  chitosan. 


Figure  9.  Contact  angles  of  2%  chitosan/0.2%  acetic  acid  and  (a)  25%  PEG  (109°) 
and  (b)  50%  PEG  (104°). 
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3.2.2  1%  Chitosan/2.0%  Acetic  Acid  Samples 


As  received,  some  of  the  amines  on  the  chitosan  are  acetylated,  resulting  in  decreased  solubility 
and  reactivity.  The  addition  of  the  acid  deacetylates  the  chitosan,  enhancing  the  water  solubility. 
As  with  the  functionalization,  there  is  likely  a  critical  concentration  of  acid  required  for 
solubility  that  does  not  require  a  large  degree  of  deacetylation.  It  is  possible  that  at  this  critical 
acid  concentration  the  chitosan  is  soluble  but  largely  nonreactive. 

To  determine  the  impact  of  acid  concentration,  the  chitosan  solution  concentration  was  decreased 
to  1%  and  the  acetic  acid  concentration  was  increased  to  2.0%.  The  solutions  were  reacted  with 
25%,  50%,  and  75%  PEG,  respectively.  All  three  samples  exhibited  a  decrease  in  the  contact 
angle,  indicating  PEG  was  present  in  the  samples.  The  25%  PEG  sample  had  a  contact  angle  of 
51.32°  ±  9.05°,  the  50%  PEG  sample  had  18.55°  ±  9.24°,  and  the  75%  PEG  sample  had  less  than 
10.93°  (figures  10a,  10b,  and  10c,  respectively).  The  75%  PEG  sample  had  droplets  of  water 
that  absorbed  too  quickly  onto  the  surface  of  the  film  and  could  not  be  measured  precisely.  An 
increase  in  the  hydrophilicity  trend  of  all  three  samples  suggests  that  there  is  a  linear  relationship 
between  the  concentration  of  PEG  and  the  degree  of  PEG  coupling  to  chitosan. 

I  I  I 


Figure  10.  Contact  angles  of  1%  chitosan/2.0%  acetic  acid  and  (a)  25%  PEG  (51°),  (b)  50%  PEG  (19°), 
and  (c)  75%  PEG  (11°). 

The  decrease  in  contact  angle,  relative  to  the  previous  samples,  may  be  attributed  to  the  dilution 
of  the  chitosan,  the  increased  acid  concentration,  or  a  cooperative  impact.  The  chitosan  dilution 
may  lead  to  a  reduction  in  the  steric  hindrance  of  neighboring  chitosan  molecules,  allowing  the 
PEG  to  react  more  quickly.  However,  it  is  anticipated  that  even  the  more  concentrated  solution 
would  have  exhibited  some  change  in  contact  angle  at  the  highest  PEG  concentration.  The 
potential  cooperative  impact  of  chitosan  dilution  and  increased  acid  concentration  will  be 
explored  by  performing  the  experiments  using  a  2%  chitosan  solution  in  2%  acetic  acid. 

FTIR  scans  of  the  samples  support  the  results  of  the  contact  angle  measurements  (figure  11).  The 
samples  containing  50%  and  75%  PEG  samples  show  an  increase  in  the  carbonyl  peak  and 
carbon-hydrogen  stretching  peak  with  increasing  peak  area.  However,  the  carbon-hydrogen 
stretching  peak  shape  changes  at  75%  PEG.  In  addition,  the  sample  containing  25%  PEG  does 
exhibit  any  significant  changes  from  the  pure  chitosan  scan  despite  a  major  change  in  contact 
angle.  This  data  may  indicate  the  limitations  of  FTIR  for  this  analysis.  It  is  possible  that  the  PEG 
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concentration  required  to  alter  the  contact  angle  is  below  the  detectable  limit  of  FTIR.  The 
ability  to  detect  the  PEG  in  low  concentration  may  also  be  limited  by  sample  geometry  as  the 
signal  is  dependent  on  the  film  placement  during  FTIR  and  the  contact  area  between  the  FTIR 
crystal  and  the  film. 

Furthermore,  it  is  counterintuitive  that  the  contact  angle  would  be  reduced  during  the  covalent 
addition  of  the  PEG  relative  to  the  blends  as  the  reaction  did  not  likely  result  in  a  quantitative 
yield.  It  is  possible  that  this  is  the  result  of  increased  deacetylation  of  the  chitosan  due  to  the 
increase  acid  concentration.  The  reduced  ability  of  the  low  surface  tension  groups  to  segregate 
to  the  air  interface  during  film  casting  might  also  be  another  factor.  This  concept  can  also  be 
explored  using  angle-resolved  XPS  and  synchrotron  NEXAFS  to  examine  the  surface 
chemistries. 


Figure  11.  FTIR  scans  of  1%  chitosan/2.0%  acetic  acid  and  25%  PEG  (purple),  50%  PEG  (blue), 
or  75%  PEG  (red). 


3.3  PPO  Samples 

3.3.1  2%  Chitosan/0.2%  Acetic  Acid  Samples 

To  determine  the  potential  impact  of  PPO  on  the  hydrophilicity  of  chitosan,  films  consisting  of 
chitosan-PPO  blends  were  produced  where  the  two  components  are  mixed  but  not  reacted.  All 
three  concentrations  of  PPO  exhibited  similar  decreases  in  contact  angle.  The  2%  chitosan/0.2% 
acetic  acid/25%  PPO  sample  had  a  contact  angle  of  34.25°  ±  5.44°  (figure  12a),  the  50%  PPO 
sample  had  34.08°  ±  9.59°  (figure  12b),  and  the  75%  PPO  blend  had  31.24°  ±  10.05° 

(figure  12c).  Similarly,  the  FTIR  scans  of  the  three  samples  closely  align,  as  seen  in  figure  13, 
with  an  increase  in  the  carbon-hydrogen  stretching  at  2900  nm"1,  This  indicates  that  PPO  has  the 
maximum  impact  at  25%  loading  and  there  is  no  added  benefit  to  higher  concentrations. 
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Figure  12.  Contact  angles  of  2%  chitosan/0.2%  acetic  acid  and  (a)  25%  PPO  (34°),  (b)  50%  PPO  (34°), 
and  (c)  75%  PPO  (31°). 


Figure  13.  FTIR  scans  of  2%  chitosan/0.2%  acetic  acid  blends  with  (a)  25%  PPO  and 
(b)  50%  PPO  (red)  and  75%  PPO  (purple). 
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3.3.2  1%  Chitosan/2.0%  Acetic  Acid  Samples 


Next,  1%  chitosan/2.0%  acetic  acid  samples  were  created  to  see  if  there  was  variability  in  the 
wettability  depending  on  the  concentration  of  acid  or  chitosan.  All  three  samples  had  a  decrease 
in  hydrophobicity;  as  seen  in  figure  14,  the  25%  PPO  sample  had  a  contact  angle  of  48.84°  ± 
21.67°,  while  the  50%  PPO  and  75%  PPO  samples  exhibited  a  contact  angle  less  than  10°,  but 
they  could  not  be  precisely  measured. 


Figure  14.  Contact  angles  of  1%  chitosan/2.0%  acetic  acid  and  (a)  25%  PPO  (49°),  (b)  50%  PPO  (<10°),  and 
(c)  75%  PPO  (<10°). 


The  FTIR  results  shown  in  figure  15  indicate  that  the  three  samples  had  similar  amounts  of  PPO 
in  each  sample.  This  result  may  further  support  that  there  are  some  limitations  of  FTIR  for  this 
project.  While  the  FTIR  does  not  exhibit  any  difference  between  the  samples,  the  contact  angle 
demonstrates  the  anticipated  trend.  Future  work  will  be  done  to  explore  the  viability  of  using 
nuclear  magnetic  resonance  (NMR)  to  determine  the  extent  of  functionality.  NMR  analysis  also 
has  some  significant  challenges  in  regards  to  sample  preparation — it  is  a  solution  technique 
rather  than  using  the  solid  film  samples  employed  in  FTIR. 


Figure  15.  FTIR  scans  of  1%  chitosan/2.0%  acetic  acid  and  25%  PPO  (green),  50%  PPO  (red), 
and  75%  (blue). 
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3.4  Epoxy  Dodecane  Samples 


The  initial  epoxy  dodecane  reaction  consisting  of  1%  chitosan  in  0.2%  acetic  acid  and  50% 
epoxy  dodecane  increased  the  hydrophilicity  to  84.42°  ±  6.49°  (figure  16).  This  contradicts  the 
theory  that  epoxy  dodecane  would  make  chitosan  more  hydrophobic  because  epoxy  dodecane 
contains  an  excessive  amount  of  hydrophobic  methylene  groups  in  its  backbone.  Additionally, 
the  FTIR  scan  did  not  show  any  indication  of  epoxy  dodecane  present  in  the  sample,  as  there  is 
no  change  from  its  FTIR  scan  and  that  of  pure  chitosan  (figure  17).  The  ring  opening  reaction 
for  the  epoxy  dodecane  may  not  have  occurred  if  the  energy  in  the  system  was  not  sufficient  to 
jumpstart  the  reaction.  In  addition,  the  epoxy  dodecane  is  not  soluble  in  the  acidified  water  and 
requires  70%  ethanol  to  form  a  uniform  mixture.  While  the  chitosan  does  remain  in  solution,  it 
is  unclear  what  impact  the  ethanol  has  on  the  reactivity.  As  a  result,  the  reaction  temperature  and 
time  may  have  to  be  modified.  In  addition,  these  experiments  need  to  be  reproduced  using  the 
lower  chitosan  concentration  and  increased  acid  concentrations  that  have  shown  promise  in 
previously  described  experiments. 
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Figure  16.  Contact  angle  of  1%  chitosan/0.2%  acetic 
acid  /50%  epoxy  dodecane  (84°). 


Figure  17.  FTIR  scan  of  1%  chitosan/0.2%  acetic  acid/50%  epoxy  dodecane,  which  has 
no  differences  compared  to  pure  chitosan. 
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4.  Summary/Conclusions 


At  this  point  in  the  project,  the  important  experimental  parameters  are  beginning  to  become 
apparent  but  will  require  further  study  to  be  conclusive.  From  the  contact  angle  data,  the  acid 
concentration  of  the  chitosan  solution  is  very  important.  An  acetic  acid  concentration  that  is  too 
low  will  solubilize  the  chitosan  but  appears  to  hinder  the  chitosan  reactivity.  Flowever,  an 
increased  acid  concentration  will  enhance  the  reactivity,  leading  to  significant  changes  in  the 
contact  angle.  In  many  cases,  FTIR  supports  the  presence  of  chemical  modifications  but  is  not 
always  quantitative.  This  is  most  prevalent  at  low  functional  densities  where  FTIR  does  not 
exhibit  a  chemical  change  in  contrast  to  a  large  change  in  contact  angle.  This  indicates  that  FTIR 
may  not  be  the  best  tool  to  monitor  the  coupling  reaction.  Other  methods  for  monitoring  the 
reaction  including  NMR  and  thin  layer  chromatography  (TLC)  are  currently  being  explored. 

While  the  changes  are  apparent  in  the  modified  chitosan  contact  angle,  there  are  a  few 
inconsistencies  that  require  further  experimentation.  Contrary  to  the  anticipated  outcome,  the 
PPO  exhibited  a  lower  contact  angle  than  PEO.  We  anticipate  that  this  is  due  to  the  PPO  being 
di-functional.  In  the  presence  of  water,  the  unreacted  epoxy  ring  will  open,  resulting  in  a 
hydroxyl  termination.  The  ability  of  the  hydroxyl  group  to  hydrogen  bond  with  the  water  may 
lead  to  the  decreased  contact  angle.  Depending  on  the  ring  opening  timescale,  it  is  possible  that 
the  blends  still  exhibited  epoxy  functionality  leading  to  an  increased  contact  angle  (-30°).  To 
determine  the  impact  of  the  end  group,  a  polymer  with  increased  chain  length  can  be  used  to 
minimize  the  impact  of  the  end  group.  In  addition,  the  impact  of  acid  concentration  and  chitosan 
dilution  needs  to  be  further  explored  to  identify  optimum  conditions. 

In  the  near  future,  we  will  also  explore  the  development  of  functional  scaffold  materials, 
including  epoxy  foams  and  soft  polymer  monoliths.  One  option  is  to  couple  the  chitosan  with 
glutaraldehyde,  another  common  biocompatible  material,  which  will  effectively  cross-link  the 
chitosan,  producing  a  polymer  gel.  These  materials  can  then  be  further  functionalized  to  provide 
specific  interactions  depending  on  the  application.  We  anticipate  that  these  methods  can  be 
readily  tailored,  scaled  up,  and  implemented  into  a  variety  of  Army  applications. 
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